1NOTE: THIS IS NOT YET CORRECT, BUT REPRESENTS THE STATE OF THE TEST RUN 2AS OF THE DATE OF THIS FILE. 3 4Fri Mar 6 11:46:58 PST 1998 5REDUCE Development Version, 28-Feb-98 ... 6 71: 1: 82: 2: 93: 3: 3: 3: 3: 3: 3: 3: 3: 10nil 11 124: 4: Comment This is a standard test file for REDUCE that has been used for 13many years. It only tests a limited number of facilities in the 14current system. In particular, it does not test floating point 15arithmetic, or any of the more advanced packages that have been made 16available since REDUCE 3.0 was released. It has been used for a long 17time to benchmark the performance of REDUCE. A description of this 18benchmarking with statistics for REDUCE 3.2 was reported in Jed B. 19Marti and Anthony C. Hearn, "REDUCE as a Lisp Benchmark", SIGSAM Bull. 2019 (1985) 8-16. That paper also gives information on the the parts of 21the system exercised by the test file. Updated statistics may be found 22in the "timings" file in the REDUCE Network Library; 23 24 25showtime; 26 27 28Time: 0 ms 29 30 31on reduce4; 32 33 34 35of type: noval 36 % For the time being. 37 38comment some examples of the FOR statement; 39 40 41comment summing the squares of the even positive integers 42 through 50; 43 44 45for i:=2 step 2 until 50 sum i**2; 46 47 4822100 49 50of type: bint 51 52 53comment to set w to the factorial of 10; 54 55 56w := for i:=1:10 product i; 57 58 593628800 60 61of type: bint 62 63 64comment alternatively, we could set the elements a(i) of the 65 array a to the factorial of i by the statements; 66 67 68array a(10); 69 70 71 72of type: noval 73 74 75a(0):=1$ 76 77 78of type: nzint 79 80 81for i:=1:10 do a(i):=i*a(i-1); 82 83 84nil 85 86of type: variable 87 88 89comment the above version of the FOR statement does not return 90 an algebraic value, but we can now use these array 91 elements as factorials in expressions, e. g.; 92 93 941+a(5); 95 96 97121 98 99of type: nzint 100 101 102comment we could have printed the values of each a(i) 103 as they were computed by writing the FOR statement as; 104 105 106for i:=1:10 do write "a(",i,") := ",a(i):= i*a(i-1); 107 108a(1) := 1 109 110a(2) := 2 111 112a(3) := 6 113 114a(4) := 24 115 116a(5) := 120 117 118a(6) := 720 119 120a(7) := 5040 121 122a(8) := 40320 123 124a(9) := 362880 125 126a(10) := 3628800 127 128 129nil 130 131of type: variable 132 133 134comment another way to use factorials would be to introduce an 135operator FAC by an integer procedure as follows; 136 137 138procedure fac(n:int) 139 begin local m:int; 140 m:=1; 141 l1: if n=0 then return m; 142 m:=m*n; 143 n:=n-1; 144 go to l1 145 end; 146 147 148fac 149 150of type: variable 151 152 153comment we can now use fac as an operator in expressions, e. g.; 154 155 156z**2+fac(4)-2*fac 2*y; 157 158 159 2 160 - 4*y + z + 24 161 162of type: xpoly 163 164 165comment note in the above example that the parentheses around 166the arguments of FAC may be omitted since it is a unary operator; 167 168 169comment the following examples illustrate the solution of some 170 complete problems; 171 172 173comment the f and g series (ref Sconzo, P., Leschack, A. R. and 174 Tobey, R. G., Astronomical Journal, Vol 70 (May 1965); 175 176 177deps:= -sigma*(mu+2*epsilon)$ 178 179 180of type: xpoly 181 182dmu:= -3*mu*sigma$ 183 184 185of type: xpoly 186 187dsig:= epsilon-2*sigma**2$ 188 189 190of type: xpoly 191 192f1:= 1$ 193 194 195of type: nzint 196 197g1:= 0$ 198 199 200of type: zero 201 202 203for i:= 1:8 do 204 <<f2:= -mu*g1 + deps*df(f1,epsilon) + dmu*df(f1,mu) + dsig*df(f1,sigma); 205 write "f(",i,") := ",f2; 206 g2:= f1 + deps*df(g1,epsilon) + dmu*df(g1,mu) + dsig*df(g1,sigma); 207 write "g(",i,") := ",g2; 208 f1:=f2; 209 g1:=g2>>; 210 211f(1) := 0 212 213g(1) := 1 214 215f(2) := - mu 216 217g(2) := 0 218 219f(3) := 3*mu*sigma 220 221g(3) := - mu 222 223 2 224f(4) := mu*(3*epsilon + mu - 15*sigma ) 225 226g(4) := 6*mu*sigma 227 228 2 229f(5) := 15*mu*sigma*( - 3*epsilon - mu + 7*sigma ) 230 231 2 232g(5) := mu*(9*epsilon + mu - 45*sigma ) 233 234 2 2 2 235f(6) := mu*( - 45*epsilon - 24*epsilon*mu + 630*epsilon*sigma - mu 236 237 2 4 238 + 210*mu*sigma - 945*sigma ) 239 240 2 241g(6) := 30*mu*sigma*( - 6*epsilon - mu + 14*sigma ) 242 243 2 2 2 244f(7) := 63*mu*sigma*(25*epsilon + 14*epsilon*mu - 150*epsilon*sigma + mu 245 246 2 4 247 - 50*mu*sigma + 165*sigma ) 248 249 2 2 2 250g(7) := mu*( - 225*epsilon - 54*epsilon*mu + 3150*epsilon*sigma - mu 251 252 2 4 253 + 630*mu*sigma - 4725*sigma ) 254 255 3 2 2 2 256f(8) := mu*(1575*epsilon + 1107*epsilon *mu - 42525*epsilon *sigma 257 258 2 2 4 259 + 117*epsilon*mu - 24570*epsilon*mu*sigma + 155925*epsilon*sigma 260 261 3 2 2 4 6 262 + mu - 2205*mu *sigma + 51975*mu*sigma - 135135*sigma ) 263 264 2 2 2 265g(8) := 126*mu*sigma*(75*epsilon + 24*epsilon*mu - 450*epsilon*sigma + mu 266 267 2 4 268 - 100*mu*sigma + 495*sigma ) 269 270 271nil 272 273of type: variable 274 275 276comment a problem in Fourier analysis; 277 278 279factor cos,sin; 280 281 282 283of type: noval 284 285 286on list; 287 288 289 290of type: noval 291 292 293(a1*cos(omega*t) + a3*cos(3*omega*t) + b1*sin(omega*t) + 294 b3*sin(3*omega*t))**3 where 295 {cos(~x)*cos(~y) => (cos(x+y)+cos(x-y))/2, 296 cos(~x)*sin(~y) => (sin(x+y)-sin(x-y))/2, 297 sin(~x)*sin(~y) => (cos(x-y)-cos(x+y))/2, 298 cos(~x)**2 => (1+cos(2*x))/2, 299 sin(~x)**2 => (1-cos(2*x))/2}; 300 301 302 3 303(3*cos(omega*t)*(a1 304 305 2 306 +a1 *a3 307 308 2 309 +2*a1*a3 310 311 2 312 +a1*b1 313 314 +2*a1*b1*b3 315 316 2 317 +2*a1*b3 318 319 2 320 -a3*b1 ) 321 322 2 323 +cos(9*omega*t)*a3*(a3 324 325 2 326 -3*b3 ) 327 328 2 329 +3*cos(7*omega*t)*(a1*a3 330 331 2 332 -a1*b3 333 334 -2*a3*b1*b3) 335 336 2 337 +3*cos(5*omega*t)*(a1 *a3 338 339 2 340 +a1*a3 341 342 -2*a1*b1*b3 343 344 2 345 -a1*b3 346 347 2 348 -a3*b1 349 350 +2*a3*b1*b3) 351 352 3 353 +cos(3*omega*t)*(a1 354 355 2 356 +6*a1 *a3 357 358 2 359 -3*a1*b1 360 361 3 362 +3*a3 363 364 2 365 +6*a3*b1 366 367 2 368 +3*a3*b3 ) 369 370 2 371 +3*sin(omega*t)*(a1 *b1 372 373 2 374 +a1 *b3 375 376 -2*a1*a3*b1 377 378 2 379 +2*a3 *b1 380 381 3 382 +b1 383 384 2 385 -b1 *b3 386 387 2 388 +2*b1*b3 ) 389 390 2 391 +sin(9*omega*t)*b3*(3*a3 392 393 2 394 -b3 ) 395 396 +3*sin(7*omega*t)*(2*a1*a3*b3 397 398 2 399 +a3 *b1 400 401 2 402 -b1*b3 ) 403 404 2 405 +3*sin(5*omega*t)*(a1 *b3 406 407 +2*a1*a3*b1 408 409 +2*a1*a3*b3 410 411 2 412 -a3 *b1 413 414 2 415 -b1 *b3 416 417 2 418 +b1*b3 ) 419 420 2 421 +sin(3*omega*t)*(3*a1 *b1 422 423 2 424 +6*a1 *b3 425 426 2 427 +3*a3 *b3 428 429 3 430 -b1 431 432 2 433 +6*b1 *b3 434 435 3 436 +3*b3 ))/4 437 438of type: xratpol 439 440 441remfac cos,sin; 442 443 444 445of type: noval 446 447 448off list; 449 450 451 452of type: noval 453 454 455comment end of Fourier analysis example; 456 457 458comment the following program, written in collaboration with David 459Barton and John Fitch, solves a problem in general relativity. it 460will compute the Einstein tensor from any given metric; 461 462 463on nero; 464 465 466 467of type: noval 468 469 470comment here we introduce the covariant and contravariant metrics; 471 472 473operator p1,q1,x; 474 475 476 477of type: noval 478 479 480array gg(3,3),h(3,3); 481 482 483 484of type: noval 485 486 487gg(0,0):=e**(q1(x(1)))$ 488 489 490of type: kernel 491 492gg(1,1):=-e**(p1(x(1)))$ 493 494 495of type: xpoly 496 497gg(2,2):=-x(1)**2$ 498 499 500of type: xpoly 501 502gg(3,3):=-x(1)**2*sin(x(2))**2$ 503 504 505of type: xpoly 506 507 508for i:=0:3 do h(i,i):=1/gg(i,i); 509 510 511nil 512 513of type: variable 514 515 516comment generate Christoffel symbols and store in arrays cs1 and cs2; 517 518 519array cs1(3,3,3),cs2(3,3,3); 520 521 522 523of type: noval 524 525 526for i:=0:3 do for j:=i:3 do 527 <<for k:=0:3 do 528 cs1(j,i,k) := cs1(i,j,k):=(df(gg(i,k),x(j))+df(gg(j,k),x(i)) - 529 df(gg(i,j),x(k)))/2; 530 for k:=0:3 do cs2(j,i,k):= cs2(i,j,k) := for p := 0:3 sum 531 h(k,p)*cs1(i,j,p)>>; 532 533 534nil 535 536of type: variable 537 538 539comment now compute the Riemann tensor and store in r(i,j,k,l); 540 541 542array r(3,3,3,3); 543 544 545 546of type: noval 547 548 549for i:=0:3 do for j:=i+1:3 do for k:=i:3 do 550 for l:=k+1:if k=i then j else 3 do 551 <<r(j,i,l,k) := r(i,j,k,l) := for q := 0:3 sum 552 gg(i,q)*(df(cs2(k,j,q),x(l))-df(cs2(j,l,q),x(k)) + 553 for p:=0:3 sum (cs2(p,l,q)*cs2(k,j,p) - 554 cs2(p,k,q)*cs2(l,j,p))); 555 r(i,j,l,k) := -r(i,j,k,l); 556 r(j,i,k,l) := -r(i,j,k,l); 557 if i neq k or j>l then 558 <<r(k,l,i,j) := r(l,k,j,i) := r(i,j,k,l); 559 r(l,k,i,j) := -r(i,j,k,l); 560 r(k,l,j,i) := -r(i,j,k,l)>>>>; 561 562 563nil 564 565of type: variable 566 567 568comment now compute and print the Ricci tensor; 569 570 571array ricci(3,3); 572 573 574 575of type: noval 576 577 578for i:=0:3 do for j:=0:3 do 579 write ricci(j,i) := ricci(i,j) := for p := 0:3 sum for q := 0:3 sum 580 h(p,q)*r(q,i,p,j); 581 5820 583 5840 585 5860 587 5880 589 5900 591 592 - df(p1(x(1)),x(1)) 593---------------------- 594 x(1) 595 5960 597 5980 599 6000 601 6020 603 604 p1(x(1)) 605 - df(p1(x(1)),x(1))*x(1) - 2*e + 2 606--------------------------------------------- 607 p1(x(1)) 608 2*e 609 6100 611 6120 613 6140 615 6160 617 618 2 p1(x(1)) 619 sin(x(2)) *( - df(p1(x(1)),x(1))*x(1) - 2*e + 2) 620---------------------------------------------------------- 621 p1(x(1)) 622 2*e 623 624 625nil 626 627of type: variable 628 629 630comment now compute and print the Ricci scalar; 631 632 633rs := for i:= 0:3 sum for j:= 0:3 sum h(i,j)*ricci(i,j); 634 635 636 p1(x(1)) 637 2*(df(p1(x(1)),x(1))*x(1) + e - 1) 638-------------------------------------------- 639 p1(x(1)) 2 640 e *x(1) 641 642of type: xratpol 643 644 645comment finally compute and print the Einstein tensor; 646 647 648array einstein(3,3); 649 650 651 652of type: noval 653 654 655for i:=0:3 do for j:=0:3 do 656 write einstein(i,j):=ricci(i,j)-rs*gg(i,j)/2; 657 658 q1(x(1)) p1(x(1)) 659 e *( - df(p1(x(1)),x(1))*x(1) - e + 1) 660------------------------------------------------------- 661 p1(x(1)) 2 662 e *x(1) 663 6640 665 6660 667 6680 669 6700 671 672 p1(x(1)) 673 e - 1 674--------------- 675 2 676 x(1) 677 6780 679 6800 681 6820 683 6840 685 686 df(p1(x(1)),x(1))*x(1) 687------------------------ 688 p1(x(1)) 689 2*e 690 6910 692 6930 694 6950 696 6970 698 699 2 700 df(p1(x(1)),x(1))*sin(x(2)) *x(1) 701----------------------------------- 702 p1(x(1)) 703 2*e 704 705 706nil 707 708of type: variable 709 710 711comment end of Einstein tensor program; 712 713 714clear gg,h,cs1,cs2,r,ricci,einstein; 715 716 717 718of type: noval 719 720 721comment an example using the matrix facility; 722 723 724matrix xx,yy,zz; 725 726 727***** matrix not defined for type list 728 729 730let xx= mat((a11,a12),(a21,a22)), 731 yy= mat((y1),(y2)); 732 733 734*** Please use => instead of = in rules 735 736*** Please use => instead of = in rules 737 738 739of type: noval 740 741 7422*det xx - 3*w; 743 744 745***** det not defined for type matrix 746 747 748zz:= xx**(-1)*yy; 749 750 751***** expt not defined for types matrix nzint 752 753 7541/xx**2; 755 756 757***** expt not defined for types matrix nzint 758 759 760comment end of matrix examples; 761 762 763comment a physics example; 764 765 766on div; 767 768 769 770of type: noval 771 comment this gives us output in same form as Bjorken and Drell; 772 773 774mass ki= 0, kf= 0, p1= m, pf= m; 775 776 777***** equal not defined for types variable zero 778 779 780vector eei,ef; 781 782 783***** ~vector not defined for type list 784 785 786mshell ki,kf,p1,pf; 787 788 789***** mshell not defined for type list 790 791 792let p1.eei= 0, p1.ef= 0, p1.pf= m**2+ki.kf, p1.ki= m*k,p1.kf= 793 m*kp, pf.eei= -kf.eei, pf.ef= ki.ef, pf.ki= m*kp, pf.kf= 794 m*k, ki.eei= 0, ki.kf= m*(k-kp), kf.ef= 0, eei.eei= -1, ef.ef= 795 -1; 796 797 798*** Please use => instead of = in rules 799 800*** cons declared operator 801 802***** ef ef invalid as list 803 804 805operator gp; 806 807 808 809of type: noval 810 811 812for all p let gp(p)= g(l,p)+m; 813 814 815*** Please use => instead of = in rules 816 817 818of type: noval 819 820 821comment this is just to save us a lot of writing; 822 823 824gp(pf)*(g(l,ef,eei,ki)/(2*ki.p1) + g(l,eei,ef,kf)/(2*kf.p1)) 825 826 827 828***** g not defined for types variable variable variable variable 829 * gp(p1)*(g(l,ki,eei,ef)/(2*ki.p1) + g(l,kf,ef,eei)/(2*kf.p1))$ 830 831 832***** g not defined for types variable variable variable variable 833 834 835write "The Compton cross-section is ",ws; 836 837 -2 838The Compton cross-section is 2*x(1) 839 840 - p1(x(1)) - p1(x(1)) 841*(e *df(p1(x(1)),x(1))*x(1) - e + 1) 842 843 844 845of type: noval 846 847 848comment end of first physics example; 849 850 851off div; 852 853 854 855of type: noval 856 857 858comment another physics example; 859 860 861index ix,iy,iz; 862 863 864***** index not defined for type list 865 866 867mass p1=mm,p2=mm,p3= mm,p4= mm,k1=0; 868 869 870***** equal not defined for types variable variable 871 872 873mshell p1,p2,p3,p4,k1; 874 875 876***** mshell not defined for type list 877 878 879vector qi,q2; 880 881 882***** ~vector not defined for type list 883 884 885factor mm,p1.p3; 886 887 888***** cons not defined for types variable variable 889 890 891order mm; 892 893 894***** order not defined for type list 895 896 897operator ga,gb; 898 899 900 901of type: noval 902 903 904for all p let ga(p)=g(la,p)+mm, gb(p)= g(lb,p)+mm; 905 906 907*** Please use => instead of = in rules 908 909*** Please use => instead of = in rules 910 911 912of type: noval 913 914 915ga(-p2)*g(la,ix)*ga(-p4)*g(la,iy)* (gb(p3)*g(lb,ix)*gb(qi) 916 917 918ga(-p2)*g(la,ix)*ga(-p4)*g(la,iy)*(gb(p3)*g(lb,ix)*gb(qi) $$$ 919 920 921***** Too few right parentheses 922 923 *g(lb,iz)*gb(p1)*g(lb,iy)*gb(q2)*g(lb,iz) + gb(p3) 924 925 926 927***** g not defined for types variable variable 928 *g(lb,iz)*gb(q2)*g(lb,ix)*gb(p1)*g(lb,iz)*gb(qi)*g(lb,iy))$ 929 930$*g(lb,iz)*gb(q2)*g(lb,ix)*gb(p1)*g(lb,iz)*gb(qi)*g(lb,iy)$$$) 931 932***** Too many right parentheses 933 934 935 936let qi=p1-k1, q2=p3+k1; 937 938 939*** Please use => instead of = in rules 940 941*** Please use => instead of = in rules 942 943 944of type: noval 945 946 947comment it is usually faster to make such substitutions after all the 948 trace algebra is done; 949 950 951write "CXN =",ws; 952 953 p1(x(1)) 954 2*(df(p1(x(1)),x(1))*x(1) + e - 1) 955CXN =-------------------------------------------- 956 p1(x(1)) 2 957 e *x(1) 958 959 960 961of type: noval 962 963 964comment end of second physics example; 965 966 967showtime; 968 969 970Time: 2346 ms plus GC time: 51 ms 971 972 973of type: noval 974 975 976end; 977 978 979of type: noval 980 9815: 5: 5: 5: 5: 5: 5: 982if(*_yyy_*:=gctime()+-*_yyy_*)>0 $$$ 983 984 985***** ; invalid in if statement 986 9876: 988;, plus GC time: prin2<<$$$then 989 990***** Improper delimiter 991 992 993*_yyy_**_yyy_* 994 995of type: variable 996 9978: 998 ms ms 999 1000of type: string 1001 1002 1003$>>terpri()>$$$> 1004 1005***** Improper delimiter 1006 100710: 10: 1008Quitting 1009Fri Mar 6 11:47:01 PST 1998 1010