Artifact 28f902f860cab1ccaaefa842f165f7637871775c3bca798436463c4d3298b5ef:
- Executable file
r38/log/alg.rlg
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 15856) [annotate] [blame] [check-ins using] [more...]
Tue Apr 15 00:32:20 2008 run on win32 Comment This is a standard test file for REDUCE that has been used for many years. It only tests a limited number of facilities in the current system. In particular, it does not test floating point arithmetic, or any of the more advanced packages that have been made available since REDUCE 3.0 was released. It does however test more than just the alg package in which it is now stored. It has been used for a long time to benchmark the performance of REDUCE. A description of this benchmarking with statistics for REDUCE 3.2 was reported in Jed B. Marti and Anthony C. Hearn, "REDUCE as a Lisp Benchmark", SIGSAM Bull. 19 (1985) 8-16. That paper also gives information on the the parts of the system exercised by the test file. Updated statistics may be found in the "timings" file in the REDUCE Network Library; showtime; Time: 0 ms comment some examples of the FOR statement; comment summing the squares of the even positive integers through 50; for i:=2 step 2 until 50 sum i**2; 22100 comment to set w to the factorial of 10; w := for i:=1:10 product i; w := 3628800 comment alternatively, we could set the elements a(i) of the array a to the factorial of i by the statements; array a(10); a(0):=1$ for i:=1:10 do a(i):=i*a(i-1); comment the above version of the FOR statement does not return an algebraic value, but we can now use these array elements as factorials in expressions, e. g.; 1+a(5); 121 comment we could have printed the values of each a(i) as they were computed by writing the FOR statement as; for i:=1:10 do write a(i):= i*a(i-1); a(1) := 1 a(2) := 2 a(3) := 6 a(4) := 24 a(5) := 120 a(6) := 720 a(7) := 5040 a(8) := 40320 a(9) := 362880 a(10) := 3628800 comment another way to use factorials would be to introduce an operator FAC by an integer procedure as follows; integer procedure fac (n); begin integer m; m:=1; l1: if n=0 then return m; m:=m*n; n:=n-1; go to l1 end; fac comment we can now use fac as an operator in expressions, e. g.; z**2+fac(4)-2*fac 2*y; 2 - 4*y + z + 24 comment note in the above example that the parentheses around the arguments of FAC may be omitted since it is a unary operator; comment the following examples illustrate the solution of some complete problems; comment the f and g series (ref Sconzo, P., Leschack, A. R. and Tobey, R. G., Astronomical Journal, Vol 70 (May 1965); deps:= -sigma*(mu+2*epsilon)$ dmu:= -3*mu*sigma$ dsig:= epsilon-2*sigma**2$ f1:= 1$ g1:= 0$ for i:= 1:8 do <<f2:=-mu*g1 + deps*df(f1,epsilon) + dmu*df(f1,mu) + dsig*df(f1,sigma); write "F(",i,") := ",f2; g2:= f1 + deps*df(g1,epsilon) + dmu*df(g1,mu) + dsig*df(g1,sigma); write "G(",i,") := ",g2; f1:=f2; g1:=g2>>; F(1) := 0 G(1) := 1 F(2) := - mu G(2) := 0 F(3) := 3*mu*sigma G(3) := - mu 2 F(4) := mu*(3*epsilon + mu - 15*sigma ) G(4) := 6*mu*sigma 2 F(5) := 15*mu*sigma*( - 3*epsilon - mu + 7*sigma ) 2 G(5) := mu*(9*epsilon + mu - 45*sigma ) 2 2 2 F(6) := mu*( - 45*epsilon - 24*epsilon*mu + 630*epsilon*sigma - mu 2 4 + 210*mu*sigma - 945*sigma ) 2 G(6) := 30*mu*sigma*( - 6*epsilon - mu + 14*sigma ) 2 2 2 F(7) := 63*mu*sigma*(25*epsilon + 14*epsilon*mu - 150*epsilon*sigma + mu 2 4 - 50*mu*sigma + 165*sigma ) 2 2 2 G(7) := mu*( - 225*epsilon - 54*epsilon*mu + 3150*epsilon*sigma - mu 2 4 + 630*mu*sigma - 4725*sigma ) 3 2 2 2 F(8) := mu*(1575*epsilon + 1107*epsilon *mu - 42525*epsilon *sigma 2 2 4 + 117*epsilon*mu - 24570*epsilon*mu*sigma + 155925*epsilon*sigma 3 2 2 4 6 + mu - 2205*mu *sigma + 51975*mu*sigma - 135135*sigma ) 2 2 2 G(8) := 126*mu*sigma*(75*epsilon + 24*epsilon*mu - 450*epsilon*sigma + mu 2 4 - 100*mu*sigma + 495*sigma ) comment a problem in Fourier analysis; factor cos,sin; on list; (a1*cos(omega*t) + a3*cos(3*omega*t) + b1*sin(omega*t) + b3*sin(3*omega*t))**3 where {cos(~x)*cos(~y) => (cos(x+y)+cos(x-y))/2, cos(~x)*sin(~y) => (sin(x+y)-sin(x-y))/2, sin(~x)*sin(~y) => (cos(x-y)-cos(x+y))/2, cos(~x)**2 => (1+cos(2*x))/2, sin(~x)**2 => (1-cos(2*x))/2}; 3 (3*cos(omega*t)*(a1 2 +a1 *a3 2 +2*a1*a3 2 +a1*b1 +2*a1*b1*b3 2 +2*a1*b3 2 -a3*b1 ) 2 +cos(9*omega*t)*a3*(a3 2 -3*b3 ) 2 +3*cos(7*omega*t)*(a1*a3 2 -a1*b3 -2*a3*b1*b3) 2 +3*cos(5*omega*t)*(a1 *a3 2 +a1*a3 -2*a1*b1*b3 2 -a1*b3 2 -a3*b1 +2*a3*b1*b3) 3 +cos(3*omega*t)*(a1 2 +6*a1 *a3 2 -3*a1*b1 3 +3*a3 2 +6*a3*b1 2 +3*a3*b3 ) 2 +3*sin(omega*t)*(a1 *b1 2 +a1 *b3 -2*a1*a3*b1 2 +2*a3 *b1 3 +b1 2 -b1 *b3 2 +2*b1*b3 ) 2 +sin(9*omega*t)*b3*(3*a3 2 -b3 ) +3*sin(7*omega*t)*(2*a1*a3*b3 2 +a3 *b1 2 -b1*b3 ) 2 +3*sin(5*omega*t)*(a1 *b3 +2*a1*a3*b1 +2*a1*a3*b3 2 -a3 *b1 2 -b1 *b3 2 +b1*b3 ) 2 +sin(3*omega*t)*(3*a1 *b1 2 +6*a1 *b3 2 +3*a3 *b3 3 -b1 2 +6*b1 *b3 3 +3*b3 ))/4 remfac cos,sin; off list; comment end of Fourier analysis example; comment the following program, written in collaboration with David Barton and John Fitch, solves a problem in general relativity. it will compute the Einstein tensor from any given metric; on nero; comment here we introduce the covariant and contravariant metrics; operator p1,q1,x; array gg(3,3),h(3,3); gg(0,0):=e**(q1(x(1)))$ gg(1,1):=-e**(p1(x(1)))$ gg(2,2):=-x(1)**2$ gg(3,3):=-x(1)**2*sin(x(2))**2$ for i:=0:3 do h(i,i):=1/gg(i,i); comment generate Christoffel symbols and store in arrays cs1 and cs2; array cs1(3,3,3),cs2(3,3,3); for i:=0:3 do for j:=i:3 do <<for k:=0:3 do cs1(j,i,k) := cs1(i,j,k):=(df(gg(i,k),x(j))+df(gg(j,k),x(i)) -df(gg(i,j),x(k)))/2; for k:=0:3 do cs2(j,i,k):= cs2(i,j,k) := for p := 0:3 sum h(k,p)*cs1(i,j,p)>>; comment now compute the Riemann tensor and store in r(i,j,k,l); array r(3,3,3,3); for i:=0:3 do for j:=i+1:3 do for k:=i:3 do for l:=k+1:if k=i then j else 3 do <<r(j,i,l,k) := r(i,j,k,l) := for q := 0:3 sum gg(i,q)*(df(cs2(k,j,q),x(l))-df(cs2(j,l,q),x(k)) + for p:=0:3 sum (cs2(p,l,q)*cs2(k,j,p) -cs2(p,k,q)*cs2(l,j,p))); r(i,j,l,k) := -r(i,j,k,l); r(j,i,k,l) := -r(i,j,k,l); if i neq k or j>l then <<r(k,l,i,j) := r(l,k,j,i) := r(i,j,k,l); r(l,k,i,j) := -r(i,j,k,l); r(k,l,j,i) := -r(i,j,k,l)>>>>; comment now compute and print the Ricci tensor; array ricci(3,3); for i:=0:3 do for j:=0:3 do write ricci(j,i) := ricci(i,j) := for p := 0:3 sum for q := 0:3 sum h(p,q)*r(q,i,p,j); q1(x(1)) ricci(0,0) := ricci(0,0) := (e *(df(p1(x(1)),x(1))*df(q1(x(1)),x(1))*x(1) 2 - 2*df(q1(x(1)),x(1),2)*x(1) - df(q1(x(1)),x(1)) *x(1) p1(x(1)) - 4*df(q1(x(1)),x(1))))/(4*e *x(1)) ricci(1,1) := ricci(1,1) := ( - df(p1(x(1)),x(1))*df(q1(x(1)),x(1))*x(1) 2 - 4*df(p1(x(1)),x(1)) + 2*df(q1(x(1)),x(1),2)*x(1) + df(q1(x(1)),x(1)) *x(1) )/(4*x(1)) ricci(2,2) := ricci(2,2) := p1(x(1)) - df(p1(x(1)),x(1))*x(1) + df(q1(x(1)),x(1))*x(1) - 2*e + 2 ---------------------------------------------------------------------- p1(x(1)) 2*e 2 ricci(3,3) := ricci(3,3) := (sin(x(2)) p1(x(1)) *( - df(p1(x(1)),x(1))*x(1) + df(q1(x(1)),x(1))*x(1) - 2*e + 2))/(2 p1(x(1)) *e ) comment now compute and print the Ricci scalar; rs := for i:= 0:3 sum for j:= 0:3 sum h(i,j)*ricci(i,j); 2 rs := (df(p1(x(1)),x(1))*df(q1(x(1)),x(1))*x(1) + 4*df(p1(x(1)),x(1))*x(1) 2 2 2 - 2*df(q1(x(1)),x(1),2)*x(1) - df(q1(x(1)),x(1)) *x(1) p1(x(1)) p1(x(1)) 2 - 4*df(q1(x(1)),x(1))*x(1) + 4*e - 4)/(2*e *x(1) ) comment finally compute and print the Einstein tensor; array einstein(3,3); for i:=0:3 do for j:=0:3 do write einstein(i,j):=ricci(i,j)-rs*gg(i,j)/2; q1(x(1)) p1(x(1)) e *( - df(p1(x(1)),x(1))*x(1) - e + 1) einstein(0,0) := ------------------------------------------------------- p1(x(1)) 2 e *x(1) p1(x(1)) - df(q1(x(1)),x(1))*x(1) + e - 1 einstein(1,1) := ------------------------------------------- 2 x(1) einstein(2,2) := (x(1)*(df(p1(x(1)),x(1))*df(q1(x(1)),x(1))*x(1) + 2*df(p1(x(1)),x(1)) - 2*df(q1(x(1)),x(1),2)*x(1) 2 - df(q1(x(1)),x(1)) *x(1) - 2*df(q1(x(1)),x(1))))/(4 p1(x(1)) *e ) 2 einstein(3,3) := (sin(x(2)) *x(1)*(df(p1(x(1)),x(1))*df(q1(x(1)),x(1))*x(1) + 2*df(p1(x(1)),x(1)) - 2*df(q1(x(1)),x(1),2)*x(1) 2 - df(q1(x(1)),x(1)) *x(1) - 2*df(q1(x(1)),x(1))))/(4 p1(x(1)) *e ) comment end of Einstein tensor program; clear gg,h,cs1,cs2,r,ricci,einstein; comment an example using the matrix facility; matrix xx,yy,zz; let xx= mat((a11,a12),(a21,a22)), yy= mat((y1),(y2)); 2*det xx - 3*w; 2*(a11*a22 - a12*a21 - 5443200) zz:= xx**(-1)*yy; [ - a12*y2 + a22*y1 ] [--------------------] [ a11*a22 - a12*a21 ] zz := [ ] [ a11*y2 - a21*y1 ] [------------------- ] [ a11*a22 - a12*a21 ] 1/xx**2; 2 a12*a21 + a22 mat((-------------------------------------------, 2 2 2 2 a11 *a22 - 2*a11*a12*a21*a22 + a12 *a21 - a12*(a11 + a22) -------------------------------------------), 2 2 2 2 a11 *a22 - 2*a11*a12*a21*a22 + a12 *a21 - a21*(a11 + a22) (-------------------------------------------, 2 2 2 2 a11 *a22 - 2*a11*a12*a21*a22 + a12 *a21 2 a11 + a12*a21 -------------------------------------------)) 2 2 2 2 a11 *a22 - 2*a11*a12*a21*a22 + a12 *a21 comment end of matrix examples; comment a physics example; on div; comment this gives us output in same form as Bjorken and Drell; mass ki= 0, kf= 0, p1= m, pf= m; vector eei,ef; mshell ki,kf,p1,pf; let p1.eei= 0, p1.ef= 0, p1.pf= m**2+ki.kf, p1.ki= m*k,p1.kf= m*kp, pf.eei= -kf.eei, pf.ef= ki.ef, pf.ki= m*kp, pf.kf= m*k, ki.eei= 0, ki.kf= m*(k-kp), kf.ef= 0, eei.eei= -1, ef.ef= -1; operator gp; for all p let gp(p)= g(l,p)+m; comment this is just to save us a lot of writing; gp(pf)*(g(l,ef,eei,ki)/(2*ki.p1) + g(l,eei,ef,kf)/(2*kf.p1)) * gp(p1)*(g(l,ki,eei,ef)/(2*ki.p1) + g(l,kf,ef,eei)/(2*kf.p1))$ write "The Compton cross-section is ",ws; 2 1 -1 1 -1 The Compton cross-section is 2*eei.ef + ---*k*kp + ---*k *kp - 1 2 2 comment end of first physics example; off div; comment another physics example; index ix,iy,iz; mass p1=mm,p2=mm,p3= mm,p4= mm,k1=0; mshell p1,p2,p3,p4,k1; vector qi,q2; factor mm,p1.p3; order mm; operator ga,gb; for all p let ga(p)=g(la,p)+mm, gb(p)= g(lb,p)+mm; ga(-p2)*g(la,ix)*ga(-p4)*g(la,iy)* (gb(p3)*g(lb,ix)*gb(qi)* g(lb,iz)*gb(p1)*g(lb,iy)*gb(q2)*g(lb,iz) + gb(p3)* g(lb,iz)*gb(q2)*g(lb,ix)*gb(p1)*g(lb,iz)*gb(qi)*g(lb,iy))$ let qi=p1-k1, q2=p3+k1; comment it is usually faster to make such substitutions after all the trace algebra is done; write "CXN =",ws; 4 4 2 2 CXN =32*mm *p1.p3 + 8*mm *(k1.p1 - k1.p3) - 16*mm *p1.p3 2 + 16*mm *p1.p3*( - k1.p1 + k1.p3 - p2.p4) 2 + 8*mm *( - k1.p1*p2.p4 - 2*k1.p2*k1.p4 + k1.p3*p2.p4) + 8*p1.p3*( k1.p2*p1.p4 - k1.p2*p3.p4 + k1.p4*p1.p2 - k1.p4*p2.p3 + 2*p1.p2*p3.p4 + 2*p1.p4*p2.p3) + 8*(k1.p1*p1.p2*p3.p4 + k1.p1*p1.p4*p2.p3 + 2*k1.p1*p2.p3*p3.p4 - 2*k1.p3*p1.p2*p1.p4 - k1.p3*p1.p2*p3.p4 - k1.p3*p1.p4*p2.p3) comment end of second physics example; showtime; Time: 64 ms plus GC time: 3 ms end; Time for test: 65 ms, plus GC time: 3 ms