reduce-historical/r37/packages/alg/alg.rlg

Sun Jan  3 23:46:01 MET 1999
REDUCE 3.7, 15-Jan-99 ...

1: 1: 
2: 2: 2: 2: 2: 2: 2: 2: 2: 
3: 3: 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: 430 ms


end;

4: 4: 4: 4: 4: 4: 4: 4: 4: 

Time for test: 430 ms
5: 5: 
Quitting
Sun Jan  3 23:46:04 MET 1999