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 does however test more 6than just the alg package in which it is now stored. It has been used 7for a long time to benchmark the performance of REDUCE. A description 8of this benchmarking with statistics for REDUCE 3.2 was reported in Jed 9B. Marti and Anthony C. Hearn, "REDUCE as a Lisp Benchmark", SIGSAM 10Bull. 19 (1985) 8-16. That paper also gives information on the the 11parts of the system exercised by the test file. Updated statistics may 12be found in the "timings" file in the REDUCE Network Library; 13 14 15showtime; 16 17 18 19 20comment some examples of the FOR statement; 21 22 23comment summing the squares of the even positive integers 24 through 50; 25 26 27for i:=2 step 2 until 50 sum i**2; 28 29 3022100 31 32 33comment to set w to the factorial of 10; 34 35 36w := for i:=1:10 product i; 37 38 39w := 3628800 40 41 42comment alternatively, we could set the elements a(i) of the 43 array a to the factorial of i by the statements; 44 45 46array a(10); 47 48 49 50a(0):=1$ 51 52 53 54for i:=1:10 do a(i):=i*a(i-1); 55 56 57 58comment the above version of the FOR statement does not return 59 an algebraic value, but we can now use these array 60 elements as factorials in expressions, e. g.; 61 62 631+a(5); 64 65 66121 67 68 69comment we could have printed the values of each a(i) 70 as they were computed by writing the FOR statement as; 71 72 73for i:=1:10 do write a(i):= i*a(i-1); 74 75 76a(1) := 1 77 78a(2) := 2 79 80a(3) := 6 81 82a(4) := 24 83 84a(5) := 120 85 86a(6) := 720 87 88a(7) := 5040 89 90a(8) := 40320 91 92a(9) := 362880 93 94a(10) := 3628800 95 96 97comment another way to use factorials would be to introduce an 98operator FAC by an integer procedure as follows; 99 100 101integer procedure fac (n); 102 begin integer m; 103 m:=1; 104 l1: if n=0 then return m; 105 m:=m*n; 106 n:=n-1; 107 go to l1 108 end; 109 110 111fac 112 113 114comment we can now use fac as an operator in expressions, e. g.; 115 116 117z**2+fac(4)-2*fac 2*y; 118 119 120 2 121 - 4*y + z + 24 122 123 124comment note in the above example that the parentheses around 125the arguments of FAC may be omitted since it is a unary operator; 126 127 128comment the following examples illustrate the solution of some 129 complete problems; 130 131 132comment the f and g series (ref Sconzo, P., Leschack, A. R. and 133 Tobey, R. G., Astronomical Journal, Vol 70 (May 1965); 134 135 136deps:= -sigma*(mu+2*epsilon)$ 137 138 139dmu:= -3*mu*sigma$ 140 141 142dsig:= epsilon-2*sigma**2$ 143 144 145f1:= 1$ 146 147 148g1:= 0$ 149 150 151 152for i:= 1:8 do 153 <<f2:=-mu*g1 + deps*df(f1,epsilon) + dmu*df(f1,mu) + dsig*df(f1,sigma); 154 write "F(",i,") := ",f2; 155 g2:= f1 + deps*df(g1,epsilon) + dmu*df(g1,mu) + dsig*df(g1,sigma); 156 write "G(",i,") := ",g2; 157 f1:=f2; 158 g1:=g2>>; 159 160 161F(1) := 0 162 163G(1) := 1 164 165F(2) := - mu 166 167G(2) := 0 168 169F(3) := 3*mu*sigma 170 171G(3) := - mu 172 173 2 174F(4) := mu*(3*epsilon + mu - 15*sigma ) 175 176G(4) := 6*mu*sigma 177 178 2 179F(5) := 15*mu*sigma*( - 3*epsilon - mu + 7*sigma ) 180 181 2 182G(5) := mu*(9*epsilon + mu - 45*sigma ) 183 184 2 2 2 185F(6) := mu*( - 45*epsilon - 24*epsilon*mu + 630*epsilon*sigma - mu 186 187 2 4 188 + 210*mu*sigma - 945*sigma ) 189 190 2 191G(6) := 30*mu*sigma*( - 6*epsilon - mu + 14*sigma ) 192 193 2 2 2 194F(7) := 63*mu*sigma*(25*epsilon + 14*epsilon*mu - 150*epsilon*sigma + mu 195 196 2 4 197 - 50*mu*sigma + 165*sigma ) 198 199 2 2 2 200G(7) := mu*( - 225*epsilon - 54*epsilon*mu + 3150*epsilon*sigma - mu 201 202 2 4 203 + 630*mu*sigma - 4725*sigma ) 204 205 3 2 2 2 206F(8) := mu*(1575*epsilon + 1107*epsilon *mu - 42525*epsilon *sigma 207 208 2 2 4 209 + 117*epsilon*mu - 24570*epsilon*mu*sigma + 155925*epsilon*sigma 210 211 3 2 2 4 6 212 + mu - 2205*mu *sigma + 51975*mu*sigma - 135135*sigma ) 213 214 2 2 2 215G(8) := 126*mu*sigma*(75*epsilon + 24*epsilon*mu - 450*epsilon*sigma + mu 216 217 2 4 218 - 100*mu*sigma + 495*sigma ) 219 220 221comment a problem in Fourier analysis; 222 223 224factor cos,sin; 225 226 227 228on list; 229 230 231 232(a1*cos(omega*t) + a3*cos(3*omega*t) + b1*sin(omega*t) 233 + b3*sin(3*omega*t))**3 234 where {cos(~x)*cos(~y) => (cos(x+y)+cos(x-y))/2, 235 cos(~x)*sin(~y) => (sin(x+y)-sin(x-y))/2, 236 sin(~x)*sin(~y) => (cos(x-y)-cos(x+y))/2, 237 cos(~x)**2 => (1+cos(2*x))/2, 238 sin(~x)**2 => (1-cos(2*x))/2}; 239 240 241 3 242(3*cos(omega*t)*(a1 243 244 2 245 +a1 *a3 246 247 2 248 +2*a1*a3 249 250 2 251 +a1*b1 252 253 +2*a1*b1*b3 254 255 2 256 +2*a1*b3 257 258 2 259 -a3*b1 ) 260 261 2 262 +cos(9*omega*t)*a3*(a3 263 264 2 265 -3*b3 ) 266 267 2 268 +3*cos(7*omega*t)*(a1*a3 269 270 2 271 -a1*b3 272 273 -2*a3*b1*b3) 274 275 2 276 +3*cos(5*omega*t)*(a1 *a3 277 278 2 279 +a1*a3 280 281 -2*a1*b1*b3 282 283 2 284 -a1*b3 285 286 2 287 -a3*b1 288 289 +2*a3*b1*b3) 290 291 3 292 +cos(3*omega*t)*(a1 293 294 2 295 +6*a1 *a3 296 297 2 298 -3*a1*b1 299 300 3 301 +3*a3 302 303 2 304 +6*a3*b1 305 306 2 307 +3*a3*b3 ) 308 309 2 310 +3*sin(omega*t)*(a1 *b1 311 312 2 313 +a1 *b3 314 315 -2*a1*a3*b1 316 317 2 318 +2*a3 *b1 319 320 3 321 +b1 322 323 2 324 -b1 *b3 325 326 2 327 +2*b1*b3 ) 328 329 2 330 +sin(9*omega*t)*b3*(3*a3 331 332 2 333 -b3 ) 334 335 +3*sin(7*omega*t)*(2*a1*a3*b3 336 337 2 338 +a3 *b1 339 340 2 341 -b1*b3 ) 342 343 2 344 +3*sin(5*omega*t)*(a1 *b3 345 346 +2*a1*a3*b1 347 348 +2*a1*a3*b3 349 350 2 351 -a3 *b1 352 353 2 354 -b1 *b3 355 356 2 357 +b1*b3 ) 358 359 2 360 +sin(3*omega*t)*(3*a1 *b1 361 362 2 363 +6*a1 *b3 364 365 2 366 +3*a3 *b3 367 368 3 369 -b1 370 371 2 372 +6*b1 *b3 373 374 3 375 +3*b3 ))/4 376 377 378remfac cos,sin; 379 380 381 382off list; 383 384 385 386comment end of Fourier analysis example; 387 388 389comment the following program, written in collaboration with David 390Barton and John Fitch, solves a problem in general relativity. it 391will compute the Einstein tensor from any given metric; 392 393 394on nero; 395 396 397 398comment here we introduce the covariant and contravariant metrics; 399 400 401operator p1,q1,x; 402 403 404 405array gg(3,3),h(3,3); 406 407 408 409gg(0,0):=e**(q1(x(1)))$ 410 411 412gg(1,1):=-e**(p1(x(1)))$ 413 414 415gg(2,2):=-x(1)**2$ 416 417 418gg(3,3):=-x(1)**2*sin(x(2))**2$ 419 420 421 422for i:=0:3 do h(i,i):=1/gg(i,i); 423 424 425 426comment generate Christoffel symbols and store in arrays cs1 and cs2; 427 428 429array cs1(3,3,3),cs2(3,3,3); 430 431 432 433for i:=0:3 do for j:=i:3 do 434 <<for k:=0:3 do 435 cs1(j,i,k) := cs1(i,j,k):=(df(gg(i,k),x(j))+df(gg(j,k),x(i)) 436 -df(gg(i,j),x(k)))/2; 437 for k:=0:3 do cs2(j,i,k):= cs2(i,j,k) := for p := 0:3 438 sum h(k,p)*cs1(i,j,p)>>; 439 440 441 442comment now compute the Riemann tensor and store in r(i,j,k,l); 443 444 445array r(3,3,3,3); 446 447 448 449for i:=0:3 do for j:=i+1:3 do for k:=i:3 do 450 for l:=k+1:if k=i then j else 3 do 451 <<r(j,i,l,k) := r(i,j,k,l) := for q := 0:3 452 sum gg(i,q)*(df(cs2(k,j,q),x(l))-df(cs2(j,l,q),x(k)) 453 + for p:=0:3 sum (cs2(p,l,q)*cs2(k,j,p) 454 -cs2(p,k,q)*cs2(l,j,p))); 455 r(i,j,l,k) := -r(i,j,k,l); 456 r(j,i,k,l) := -r(i,j,k,l); 457 if i neq k or j>l 458 then <<r(k,l,i,j) := r(l,k,j,i) := r(i,j,k,l); 459 r(l,k,i,j) := -r(i,j,k,l); 460 r(k,l,j,i) := -r(i,j,k,l)>>>>; 461 462 463 464comment now compute and print the Ricci tensor; 465 466 467array ricci(3,3); 468 469 470 471for i:=0:3 do for j:=0:3 do 472 write ricci(j,i) := ricci(i,j) := for p := 0:3 sum for q := 0:3 473 sum h(p,q)*r(q,i,p,j); 474 475 476 q1(x(1)) 477ricci(0,0) := ricci(0,0) := (e *(df(p1(x(1)),x(1))*df(q1(x(1)),x(1))*x(1) 478 479 2 480 - 2*df(q1(x(1)),x(1),2)*x(1) - df(q1(x(1)),x(1)) *x(1) 481 482 p1(x(1)) 483 - 4*df(q1(x(1)),x(1))))/(4*e *x(1)) 484 485ricci(1,1) := ricci(1,1) := ( - df(p1(x(1)),x(1))*df(q1(x(1)),x(1))*x(1) 486 487 2 488 - 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) 489 490 )/(4*x(1)) 491 492ricci(2,2) := ricci(2,2) := 493 494 p1(x(1)) 495 - df(p1(x(1)),x(1))*x(1) + df(q1(x(1)),x(1))*x(1) - 2*e + 2 496---------------------------------------------------------------------- 497 p1(x(1)) 498 2*e 499 500 2 501ricci(3,3) := ricci(3,3) := (sin(x(2)) 502 503 p1(x(1)) 504 *( - df(p1(x(1)),x(1))*x(1) + df(q1(x(1)),x(1))*x(1) - 2*e + 2))/(2 505 506 p1(x(1)) 507 *e ) 508 509 510comment now compute and print the Ricci scalar; 511 512 513rs := for i:= 0:3 sum for j:= 0:3 sum h(i,j)*ricci(i,j); 514 515 516 2 517rs := (df(p1(x(1)),x(1))*df(q1(x(1)),x(1))*x(1) + 4*df(p1(x(1)),x(1))*x(1) 518 519 2 2 2 520 - 2*df(q1(x(1)),x(1),2)*x(1) - df(q1(x(1)),x(1)) *x(1) 521 522 p1(x(1)) p1(x(1)) 2 523 - 4*df(q1(x(1)),x(1))*x(1) + 4*e - 4)/(2*e *x(1) ) 524 525 526comment finally compute and print the Einstein tensor; 527 528 529array einstein(3,3); 530 531 532 533for i:=0:3 do for j:=0:3 do 534 write einstein(i,j):=ricci(i,j)-rs*gg(i,j)/2; 535 536 537 q1(x(1)) p1(x(1)) 538 e *( - df(p1(x(1)),x(1))*x(1) - e + 1) 539einstein(0,0) := ------------------------------------------------------- 540 p1(x(1)) 2 541 e *x(1) 542 543 p1(x(1)) 544 - df(q1(x(1)),x(1))*x(1) + e - 1 545einstein(1,1) := ------------------------------------------- 546 2 547 x(1) 548 549einstein(2,2) := (x(1)*(df(p1(x(1)),x(1))*df(q1(x(1)),x(1))*x(1) 550 551 + 2*df(p1(x(1)),x(1)) - 2*df(q1(x(1)),x(1),2)*x(1) 552 553 2 554 - df(q1(x(1)),x(1)) *x(1) - 2*df(q1(x(1)),x(1))))/(4 555 556 p1(x(1)) 557 *e ) 558 559 2 560einstein(3,3) := (sin(x(2)) *x(1)*(df(p1(x(1)),x(1))*df(q1(x(1)),x(1))*x(1) 561 562 + 2*df(p1(x(1)),x(1)) - 2*df(q1(x(1)),x(1),2)*x(1) 563 564 2 565 - df(q1(x(1)),x(1)) *x(1) - 2*df(q1(x(1)),x(1))))/(4 566 567 p1(x(1)) 568 *e ) 569 570 571comment end of Einstein tensor program; 572 573 574clear gg,h,cs1,cs2,r,ricci,einstein; 575 576 577 578comment an example using the matrix facility; 579 580 581matrix xx,yy,zz; 582 583 584 585let xx= mat((a11,a12),(a21,a22)), 586 yy= mat((y1),(y2)); 587 588 589 5902*det xx - 3*w; 591 592 5932*(a11*a22 - a12*a21 - 5443200) 594 595 596zz:= xx**(-1)*yy; 597 598 599 [ - a12*y2 + a22*y1 ] 600 [--------------------] 601 [ a11*a22 - a12*a21 ] 602zz := [ ] 603 [ a11*y2 - a21*y1 ] 604 [------------------- ] 605 [ a11*a22 - a12*a21 ] 606 607 608 6091/xx**2; 610 611 612 2 613 a12*a21 + a22 614mat((-------------------------------------------, 615 2 2 2 2 616 a11 *a22 - 2*a11*a12*a21*a22 + a12 *a21 617 618 - a12*(a11 + a22) 619 -------------------------------------------), 620 2 2 2 2 621 a11 *a22 - 2*a11*a12*a21*a22 + a12 *a21 622 623 - a21*(a11 + a22) 624 (-------------------------------------------, 625 2 2 2 2 626 a11 *a22 - 2*a11*a12*a21*a22 + a12 *a21 627 628 2 629 a11 + a12*a21 630 -------------------------------------------)) 631 2 2 2 2 632 a11 *a22 - 2*a11*a12*a21*a22 + a12 *a21 633 634 635 636comment end of matrix examples; 637 638 639comment a physics example; 640 641 642on div; 643 644 comment this gives us output in same form as Bjorken and Drell; 645 646 647mass ki= 0, kf= 0, p1= m, pf= m; 648 649 650 651vector eei,ef; 652 653 654 655mshell ki,kf,p1,pf; 656 657 658 659let p1.eei= 0, p1.ef= 0, p1.pf= m**2+ki.kf, p1.ki= m*k,p1.kf= 660 m*kp, pf.eei= -kf.eei, pf.ef= ki.ef, pf.ki= m*kp, pf.kf= 661 m*k, ki.eei= 0, ki.kf= m*(k-kp), kf.ef= 0, eei.eei= -1, ef.ef= 662 -1; 663 664 665 666operator gp; 667 668 669 670for all p let gp(p)= g(l,p)+m; 671 672 673 674comment this is just to save us a lot of writing; 675 676 677gp(pf)*(g(l,ef,eei,ki)/(2*ki.p1) + g(l,eei,ef,kf)/(2*kf.p1)) 678 * gp(p1)*(g(l,ki,eei,ef)/(2*ki.p1) + g(l,kf,ef,eei)/(2*kf.p1))$ 679 680 681 682write "The Compton cross-section is ",ws; 683 684 685 2 1 -1 1 -1 686The Compton cross-section is 2*eei.ef + ---*k*kp + ---*k *kp - 1 687 2 2 688 689 690comment end of first physics example; 691 692 693off div; 694 695 696 697comment another physics example; 698 699 700index ix,iy,iz; 701 702 703 704mass p1=mm,p2=mm,p3= mm,p4= mm,k1=0; 705 706 707 708mshell p1,p2,p3,p4,k1; 709 710 711 712vector qi,q2; 713 714 715 716factor mm,p1.p3; 717 718 719 720order mm; 721 722 723 724operator ga,gb; 725 726 727 728for all p let ga(p)=g(la,p)+mm, gb(p)= g(lb,p)+mm; 729 730 731 732ga(-p2)*g(la,ix)*ga(-p4)*g(la,iy)* (gb(p3)*g(lb,ix)*gb(qi)* 733 g(lb,iz)*gb(p1)*g(lb,iy)*gb(q2)*g(lb,iz) + gb(p3)* 734 g(lb,iz)*gb(q2)*g(lb,ix)*gb(p1)*g(lb,iz)*gb(qi)*g(lb,iy))$ 735 736 737 738let qi=p1-k1, q2=p3+k1; 739 740 741 742comment it is usually faster to make such substitutions after all the 743 trace algebra is done; 744 745 746write "CXN =",ws; 747 748 749 4 4 2 2 750CXN =32*mm *p1.p3 + 8*mm *(k1.p1 - k1.p3) - 16*mm *p1.p3 751 752 2 753 + 16*mm *p1.p3*( - k1.p1 + k1.p3 - p2.p4) 754 755 2 756 + 8*mm *( - k1.p1*p2.p4 - 2*k1.p2*k1.p4 + k1.p3*p2.p4) + 8*p1.p3*( 757 758 k1.p2*p1.p4 - k1.p2*p3.p4 + k1.p4*p1.p2 - k1.p4*p2.p3 + 2*p1.p2*p3.p4 759 760 + 2*p1.p4*p2.p3) + 8*(k1.p1*p1.p2*p3.p4 + k1.p1*p1.p4*p2.p3 761 762 + 2*k1.p1*p2.p3*p3.p4 - 2*k1.p3*p1.p2*p1.p4 - k1.p3*p1.p2*p3.p4 763 764 - k1.p3*p1.p4*p2.p3) 765 766 767comment end of second physics example; 768 769 770showtime; 771 772 773 774 775end; 776 777Tested on x86_64-pc-windows CSL 778Time (counter 1): 0 ms 779 780End of Lisp run after 0.00+0.04 seconds 781real 0.22 782user 0.03 783sys 0.04 784