@@ -1,775 +1,775 @@
-REDUCE 3.6, 15-Jul-95, patched to 6 Mar 96 ...
-
-
-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 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: 20 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;
-
-
-
-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: 1030 ms  plus GC time: 10 ms
-
-end;
-(TIME:  reduce 1050 1060)
+REDUCE 3.6, 15-Jul-95, patched to 6 Mar 96 ...
+
+
+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 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: 20 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;
+
+
+
+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: 1030 ms  plus GC time: 10 ms
+
+end;
+(TIME:  reduce 1050 1060)