1COMMENT This is a standard test file for REDUCE that has been used for 2many years. It only tests a limited number of facilities in the 3current system. In particular, it does not test floating point 4arithmetic, or any of the more advanced packages that have been made 5available since REDUCE 3.0 was released. It has been used for a long 6time to benchmark the performance of REDUCE. A description of this 7benchmarking with statistics for REDUCE 3.2 was reported in Jed B. 8Marti and Anthony C. Hearn, "REDUCE as a Lisp Benchmark", SIGSAM Bull. 919 (1985) 8-16. That paper also gives information on the the parts of 10the system exercised by the test file. Updated statistics may be found 11in the "timings" file in the REDUCE Network Library; 12 13showtime; 14 15on reduce4; % For the time being. 16 17COMMENT some examples of the FOR statement; 18 19COMMENT summing the squares of the even positive integers 20 through 50; 21 22for i:=2 step 2 until 50 sum i**2; 23 24COMMENT to set w to the factorial of 10; 25 26w := for i:=1:10 product i; 27 28COMMENT alternatively, we could set the elements a(i) of the 29 array a to the factorial of i by the statements; 30 31array a(10); 32 33a(0):=1$ 34 35for i:=1:10 do a(i):=i*a(i-1); 36 37COMMENT the above version of the FOR statement does not return 38 an algebraic value, but we can now use these array 39 elements as factorials in expressions, e. g.; 40 411+a(5); 42 43COMMENT we could have printed the values of each a(i) 44 as they were computed by writing the FOR statement as; 45 46for i:=1:10 do write "a(",i,") := ",a(i):= i*a(i-1); 47 48COMMENT another way to use factorials would be to introduce an 49operator FAC by an integer procedure as follows; 50 51procedure fac(n:int) 52 begin local m:int; 53 m:=1; 54 l1: if n=0 then return m; 55 m:=m*n; 56 n:=n-1; 57 go to l1 58 end; 59 60COMMENT we can now use fac as an operator in expressions, e. g.; 61 62z**2+fac(4)-2*fac 2*y; 63 64COMMENT note in the above example that the parentheses around 65the arguments of FAC may be omitted since it is a unary operator; 66 67COMMENT the following examples illustrate the solution of some 68 complete problems; 69 70COMMENT the f and g series (ref Sconzo, P., Leschack, A. R. and 71 Tobey, R. G., Astronomical Journal, Vol 70 (May 1965); 72 73deps:= -sigma*(mu+2*epsilon)$ 74dmu:= -3*mu*sigma$ 75dsig:= epsilon-2*sigma**2$ 76f1:= 1$ 77g1:= 0$ 78 79for i:= 1:8 do 80 <<f2:= -mu*g1 + deps*df(f1,epsilon) + dmu*df(f1,mu) + dsig*df(f1,sigma); 81 write "f(",i,") := ",f2; 82 g2:= f1 + deps*df(g1,epsilon) + dmu*df(g1,mu) + dsig*df(g1,sigma); 83 write "g(",i,") := ",g2; 84 f1:=f2; 85 g1:=g2>>; 86 87COMMENT a problem in Fourier analysis; 88 89factor cos,sin; 90 91on list; 92 93(a1*cos(omega*t) + a3*cos(3*omega*t) + b1*sin(omega*t) + 94 b3*sin(3*omega*t))**3 where 95 {cos(~x)*cos(~y) => (cos(x+y)+cos(x-y))/2, 96 cos(~x)*sin(~y) => (sin(x+y)-sin(x-y))/2, 97 sin(~x)*sin(~y) => (cos(x-y)-cos(x+y))/2, 98 cos(~x)**2 => (1+cos(2*x))/2, 99 sin(~x)**2 => (1-cos(2*x))/2}; 100 101remfac cos,sin; 102 103off list; 104 105COMMENT end of Fourier analysis example; 106 107COMMENT the following program, written in collaboration with David 108Barton and John Fitch, solves a problem in general relativity. it 109will compute the Einstein tensor from any given metric; 110 111on nero; 112 113COMMENT here we introduce the covariant and contravariant metrics; 114 115operator p1,q1,x; 116 117array gg(3,3),h(3,3); 118 119gg(0,0):=e**(q1(x(1)))$ 120gg(1,1):=-e**(p1(x(1)))$ 121gg(2,2):=-x(1)**2$ 122gg(3,3):=-x(1)**2*sin(x(2))**2$ 123 124for i:=0:3 do h(i,i):=1/gg(i,i); 125 126COMMENT generate Christoffel symbols and store in arrays cs1 and cs2; 127 128array cs1(3,3,3),cs2(3,3,3); 129 130for i:=0:3 do for j:=i:3 do 131 <<for k:=0:3 do 132 cs1(j,i,k) := cs1(i,j,k):=(df(gg(i,k),x(j))+df(gg(j,k),x(i)) - 133 df(gg(i,j),x(k)))/2; 134 for k:=0:3 do cs2(j,i,k):= cs2(i,j,k) := for p := 0:3 sum 135 h(k,p)*cs1(i,j,p)>>; 136 137COMMENT now compute the Riemann tensor and store in r(i,j,k,l); 138 139array r(3,3,3,3); 140 141for i:=0:3 do for j:=i+1:3 do for k:=i:3 do 142 for l:=k+1:if k=i then j else 3 do 143 <<r(j,i,l,k) := r(i,j,k,l) := for q := 0:3 sum 144 gg(i,q)*(df(cs2(k,j,q),x(l))-df(cs2(j,l,q),x(k)) + 145 for p:=0:3 sum (cs2(p,l,q)*cs2(k,j,p) - 146 cs2(p,k,q)*cs2(l,j,p))); 147 r(i,j,l,k) := -r(i,j,k,l); 148 r(j,i,k,l) := -r(i,j,k,l); 149 if i neq k or j>l then 150 <<r(k,l,i,j) := r(l,k,j,i) := r(i,j,k,l); 151 r(l,k,i,j) := -r(i,j,k,l); 152 r(k,l,j,i) := -r(i,j,k,l)>>>>; 153 154COMMENT now compute and print the Ricci tensor; 155 156array ricci(3,3); 157 158for i:=0:3 do for j:=0:3 do 159 write ricci(j,i) := ricci(i,j) := for p := 0:3 sum for q := 0:3 sum 160 h(p,q)*r(q,i,p,j); 161 162COMMENT now compute and print the Ricci scalar; 163 164rs := for i:= 0:3 sum for j:= 0:3 sum h(i,j)*ricci(i,j); 165 166COMMENT finally compute and print the Einstein tensor; 167 168array einstein(3,3); 169 170for i:=0:3 do for j:=0:3 do 171 write einstein(i,j):=ricci(i,j)-rs*gg(i,j)/2; 172 173COMMENT end of Einstein tensor program; 174 175clear gg,h,cs1,cs2,r,ricci,einstein; 176 177COMMENT an example using the matrix facility; 178 179matrix xx,yy,zz; 180 181let xx= mat((a11,a12),(a21,a22)), 182 yy= mat((y1),(y2)); 183 1842*det xx - 3*w; 185 186zz:= xx**(-1)*yy; 187 1881/xx**2; 189 190COMMENT end of matrix examples; 191 192COMMENT a physics example; 193 194on div; COMMENT this gives us output in same form as Bjorken and Drell; 195 196mass ki= 0, kf= 0, p1= m, pf= m; 197 198vector eei,ef; 199 200mshell ki,kf,p1,pf; 201 202let p1.eei= 0, p1.ef= 0, p1.pf= m**2+ki.kf, p1.ki= m*k,p1.kf= 203 m*kp, pf.eei= -kf.eei, pf.ef= ki.ef, pf.ki= m*kp, pf.kf= 204 m*k, ki.eei= 0, ki.kf= m*(k-kp), kf.ef= 0, eei.eei= -1, ef.ef= 205 -1; 206 207operator gp; 208 209for all p let gp(p)= g(l,p)+m; 210 211COMMENT this is just to save us a lot of writing; 212 213gp(pf)*(g(l,ef,eei,ki)/(2*ki.p1) + g(l,eei,ef,kf)/(2*kf.p1)) 214 * gp(p1)*(g(l,ki,eei,ef)/(2*ki.p1) + g(l,kf,ef,eei)/(2*kf.p1))$ 215 216write "The Compton cross-section is ",ws; 217 218COMMENT end of first physics example; 219 220off div; 221 222COMMENT another physics example; 223 224index ix,iy,iz; 225 226mass p1=mm,p2=mm,p3= mm,p4= mm,k1=0; 227 228mshell p1,p2,p3,p4,k1; 229 230vector qi,q2; 231 232factor mm,p1.p3; 233 234order mm; 235 236operator ga,gb; 237 238for all p let ga(p)=g(la,p)+mm, gb(p)= g(lb,p)+mm; 239 240ga(-p2)*g(la,ix)*ga(-p4)*g(la,iy)* (gb(p3)*g(lb,ix)*gb(qi)* 241 g(lb,iz)*gb(p1)*g(lb,iy)*gb(q2)*g(lb,iz) + gb(p3)* 242 g(lb,iz)*gb(q2)*g(lb,ix)*gb(p1)*g(lb,iz)*gb(qi)*g(lb,iy))$ 243 244let qi=p1-k1, q2=p3+k1; 245 246COMMENT it is usually faster to make such substitutions after all the 247 trace algebra is done; 248 249write "CXN =",ws; 250 251COMMENT end of second physics example; 252 253showtime; 254 255end; 256