Artifact c99fe40e3615928f8e9a88b9a9862489c02bee7406fa11c8c8a905f0f128287e:
- Executable file
r38/packages/arith/arith.rlg
— part of check-in
[f2fda60abd]
at
2011-09-02 18:13:33
on branch master
— Some historical releases purely for archival purposes
git-svn-id: https://svn.code.sf.net/p/reduce-algebra/code/trunk/historical@1375 2bfe0521-f11c-4a00-b80e-6202646ff360 (user: arthurcnorman@users.sourceforge.net, size: 34575) [annotate] [blame] [check-ins using] [more...]
Tue Feb 10 12:26:19 2004 run on Linux % Tests of REDUCE Arithmetic. % Authors: Anthony C. Hearn and Stanley L. Kameny. % Copyright (c) 1987, 1988, 1989, 1991, Stanley L. Kameny. % Copyright (c) 1998, Anthony C. Hearn. % All Rights Reserved. % This test file is a combination of three tests files from versions % 3.6 or earlier: math, rounded and complex. % Simple tests of rounded arithmetic. % Tests in the exact mode. x := 1/2; 1 x := --- 2 y := x + 0.7; 6 y := --- 5 % Tests in approximate mode. on rounded; y; 1.2 % As expected not converted to approximate form. z := y+1.2; z := 2.4 z/3; 0.8 % Let's raise this to a high power. ws^24; 0.00472236648287 % Now a high exponent value. % 10.2^821; % Elementary function evaluation. cos(pi); - 1 symbolic ws; (!*sq ((!:rd!: . -1.0) . 1) t) z := sin(pi); z := 1.22464679915e-16 symbolic ws; (!*sq ((!:rd!: . 1.2246467991474e-16) . 1) t) % Handling very small quantities. % With normal defaults, underflows are converted to 0. exp(-100000.1**2); 0 % However, if you really want that small number, roundbf can be used. on roundbf; exp(-100000.1**2); 1.18441281937e-4342953505 off roundbf; % Now let us evaluate pi. pi; 3.14159265359 % Let us try a higher precision. prec0 := precision 50; prec0 := 12 pi; 3.1415926535897932384626433832795028841971693993751 % Now find the cosine of pi/6. cos(ws/6); 0.86602540378443864676372317075293618347140262690519 % This should be the sqrt(3)/2. ws**2; 0.75 % Here are some well known examples which show the power of this system. precision 10; 50 % This should give the usual default again. let xx=e**(pi*sqrt(163)); let yy=1-2*cos((6*log(2)+log(10005))/sqrt(163)); % First notice that xx looks like an integer. xx; 2.625374126e+17 % And that yy looks like zero. yy; 0 % But of course it's an illusion. precision 50; 10 xx; 2.6253741264076874399999999999925007259719818568888e+17 yy; - 1.2815256559456092775159749532170513334408547400481e-16 % Now let's look at an unusual way of finding an old friend. procedure agm; begin scalar a,b,u,x,y,p,pn; a := 1; b := 1/sqrt 2; u:= 1/4; x := 1$ pn := 4; repeat <<p := pn; y := a; a := (a+b)/2; b := sqrt(y*b); % Arith-geom mean. u := u-x*(a-y)**2; x := 2*x; pn := a**2/u; write "pn=",pn>> until pn>=p; return p end; agm agm(); pn=3.1876726427121086272019299705253692326510535718594 pn=3.1416802932976532939180704245600093827957194388154 pn=3.1415926538954464960029147588180434861088792372613 pn=3.1415926535897932384663606027066313217577024113424 pn=3.1415926535897932384626433832795028841971699491647 pn=3.1415926535897932384626433832795028841971693993751 pn=3.1415926535897932384626433832795028841971693993751 3.1415926535897932384626433832795028841971693993751 % The limit is obviously. pi; 3.1415926535897932384626433832795028841971693993751 off rounded; clear x; precision prec0; 50 % Tests of Complex arithmetic. on complex; (31+i)/74; 31 + i -------- 74 ws/(b+1); 31 + i ------------ 74*(b + 1) % This now comes out right! w:=(x+3*i)**2; 2 w := x + 6*i*x - 9 on gcd; (x**3-7*x**2+x-7)/(x**2+(3+i)*x+3*i); 2 x - (7 + i)*x + 7*i ---------------------- x + 3 off gcd; sqrt(x**4+14*i*x**3-51*x**2-14*i*x+1); 2 x + 7*i*x - 1 % All rounded tests are done twice: first, they are done at the default % precision, in which all rounded operations use standard floating point % logic. Then precision is increased, causing all rounded operations to % use extended precision bigfloat arithmetic. This is necessary to % exercise and test the bigfloat-based arithmetic functions. prec0 := precision 0; prec0 := 12 % To determine the nominal default precision. % Tests using default precision. on rounded; *** Domain mode complex changed to complex-rounded (3.25 + 8.5i) + (6.75 - 8.5i); 10.0 (3.25 + 8.5i) - (6.0 - 9.5i); - 2.75 + 18.0*i (1.0 + 10.0*i)*(-6.5 + 2.5*i); - 31.5 - 62.5*i (1.2 - 3.4*i)*(-5.6 + 7.8*i); 19.8 + 28.4*i (19.8 + 28.4*i)/(-5.6 + 7.8*i); 1.2 - 3.4*i e; 2.71828182846 pi; 3.14159265359 17*i**2; -17 (-7.0 + 24.0*i)**(1/2); 3.0 + 4.0*i sqrt(-7.0 + 24.0*i); 3.0 + 4.0*i sqrt(-10.12 - 8.16*i); 1.2 - 3.4*i sin(0.0 + 0.0*i); 0 sin(1.0 + 0.0*i); 0.841470984808 sin(1.0 + 1.0*i); 1.29845758142 + 0.634963914785*i cos(0.0 + 0.0*i); 1 cos(1.0 - 0.0*i); 0.540302305868 cos(1.0 + 1.0*i); 0.833730025131 - 0.988897705763*i tan(0.0 + 0.0*i); 0 tan(1.0 + 0.0*i); 1.55740772465 tan(1.0 + 1.0*i); 0.27175258532 + 1.08392332734*i asin(1.0 + 1.0*i); 0.666239432493 + 1.06127506191*i acos(1.0 + 1.0*i); 0.904556894302 - 1.06127506191*i atan(1.0 + 1.0*i); 1.0172219679 + 0.402359478109*i log(1.0 + 1.0*i); 0.34657359028 + 0.785398163397*i asin 2; 1.57079632679 - 1.31695789692*i sin ws; 2.0 - 1.06057523872e-16*i acos 2; 1.31695789692*i cos ws; 2.0 atan(1+i); 1.0172219679 + 0.402359478109*i tan ws; 1 + i log(2+i); 0.804718956217 + 0.463647609001*i exp ws; 2.0 + i e**(i*pi); - 1 + 1.22464679915e-16*i e**i; 0.540302305868 + 0.841470984808*i z := sqrt i; z := 0.707106781187 + 0.707106781187*i z**2; i off rounded; *** Domain mode complex-rounded changed to complex %-----------------end of normal floating point tests-------------------- precision(prec0+6); 12 % Arbitrary precision increase -> bigfloat functions %----------------------start of bigfloat tests-------------------------- on rounded; *** Domain mode complex changed to complex-rounded (3.25 + 8.5i) + (6.75 - 8.5i); 10.0 (3.25 + 8.5i) - (6.0 - 9.5i); - 2.75 + 18.0*i (1.0 + 10.0*i)*(-6.5 + 2.5*i); - 31.5 - 62.5*i (1.2 - 3.4*i)*(-5.6 + 7.8*i); 19.8 + 28.4*i (19.8 + 28.4*i)/(-5.6 + 7.8*i); 1.2 - 3.4*i e; 2.71828182845904524 pi; 3.14159265358979324 17*i**2; -17 (-7.0 + 24.0*i)**(1/2); 3.0 + 4.0*i sqrt(-7.0 + 24.0*i); 3.0 + 4.0*i sqrt(-10.12 - 8.16*i); 1.2 - 3.4*i sin(0.0 + 0.0*i); 0 sin(1.0 + 0.0*i); 0.841470984807896507 sin(1.0 + 1.0*i); 1.29845758141597729 + 0.634963914784736108*i cos(0.0 + 0.0*i); 1 cos(1.0 - 0.0*i); 0.540302305868139717 cos(1.0 + 1.0*i); 0.833730025131149049 - 0.988897705762865096*i tan(0.0 + 0.0*i); 0 tan(1.0 + 0.0*i); 1.55740772465490223 tan(1.0 + 1.0*i); 0.271752585319511717 + 1.08392332733869454*i asin(1.0 + 1.0*i); 0.666239432492515255 + 1.06127506190503565*i acos(1.0 + 1.0*i); 0.904556894302381364 - 1.06127506190503565*i atan(1.0 + 1.0*i); 1.01722196789785137 + 0.402359478108525094*i log(1.0 + 1.0*i); 0.346573590279972655 + 0.78539816339744831*i asin 2; 1.57079632679489662 - 1.31695789692481671*i sin ws; 2.0 acos 2; 1.31695789692481671*i cos ws; 2.0 atan(1+i); 1.01722196789785137 + 0.402359478108525094*i tan ws; 1 + i log(2+i); 0.804718956217050187 + 0.463647609000806116*i exp ws; 2.0 + i e**(i*pi); - 1 e**i; 0.540302305868139717 + 0.841470984807896507*i z := sqrt i; z := 0.707106781186547524 + 0.707106781186547524*i z**2; i off rounded; *** Domain mode complex-rounded changed to complex % --------------------------------------------------------------------- % The following examples are independent of precision. precision prec0; 18 % Restores default precision. s:= 1.1+2.3i; 11 + 23*i s := ----------- 10 s/4; 11 + 23*i ----------- 40 % This would have had a common factor of 4. x:= a+1.1+2.3i; 10*a + 11 + 23*i x := ------------------ 10 y:= b+1.2+1.3i; 10*b + 12 + 13*i y := ------------------ 10 z:= x/y; 100*a*b + (120 - 130*i)*a + (110 + 230*i)*b + 431 + 133*i z := ----------------------------------------------------------- 2 100*b + 240*b + 313 z/4; 100*a*b + (120 - 130*i)*a + (110 + 230*i)*b + 431 + 133*i ----------------------------------------------------------- 2 4*(100*b + 240*b + 313) % This would have had a common polynomial factor b^2 + ... z*7/4; 7*(100*a*b + (120 - 130*i)*a + (110 + 230*i)*b + 431 + 133*i) --------------------------------------------------------------- 2 4*(100*b + 240*b + 313) s/(c^2+c+1); 11 + 23*i ----------------- 2 10*(c + c + 1) % This would have had a common factor of c^2+c+1, clear x; zz:= x^2+(1.1+2.3i)*x+1.2+1.3i; 2 10*x + (11 + 23*i)*x + 12 + 13*i zz := ----------------------------------- 10 ss:=1.23456789x^2+1.3579i*x+5.6789; 2 123456789*x + 135790000*i*x + 567890000 ss := ------------------------------------------ 100000000 z:= x+1.1+2.3i; 10*x + 11 + 23*i z := ------------------ 10 on rationalize; z; 10*x + 11 + 23*i ------------------ 10 % Same as previous answer. off rationalize; 1.23456789x^2+2.3456i*x+7.89; 2 123456789*x + 234560000*i*x + 789000000 ------------------------------------------ 100000000 on factor; x**2+1; (x + i)*(x - i) x**4-1; (x + i)*(x - i)*(x + 1)*(x - 1) x**4+(i+2)*x**3+(2*i+5)*x**2+(2*i+6)*x+6; 2 (x + i*x + 3)*(x + 1 + i)*(x + 1 - i) (2*i+3)*x**4+(3*i-2)*x**3-2*(i+1)*x**2+i*x-1; 2 2 i*((2 - 3*i)*x - i)*(x + i*x - 1) % Multivariate examples. x**2+y**2; (10*b + 10*i*x + 12 + 13*i)*(10*b - 10*i*x + 12 + 13*i) --------------------------------------------------------- 100 off factor; factorize(x**2+1); {{x + i,1},{x - i,1}} off complex; % Tests of some elementary functions. comment Integer functions that work in all domain modes, independent of switch NUMVAL, so long as their arguments evaluate to real numbers. Functions of one argument: FIX, SGN, ROUND, CEILING, FLOOR (The following functions are available only in symbolic mode, so they are not tested here: ISQRT, ICBRT, ILOG2, IROOTN); fix a; fix(a) % Will be evaluated only if a evaluates to a real number. a := 27/4; 27 a := ---- 4 fix a; 6 fix 12.345; 12 sign (-15/2); -1 round 12.5; 13 ceiling 12.5; 13 floor 12.5; 12 % isqrt 12.5; % icbrt 12.5; % ilog2 130.7; % irootn(72,4); % irootn(72,3/2); % This will not evaluate. comment Functions which require arguments which evaluate to integers. Function of one argument: FACTORIAL Fumction of two arguments: PERM, CHOOSE; $ factorial 10; 3628800 perm(5,10); 30240 % Permutations of 5 out of 10. choose(5,10); 252 % Choose 5 out of 10. comment These functions are evaluated in dmodes ROUNDED and COMPLEX-ROUNDED (ON ROUNDED,COMPLEX) so long as their arguments and values evaluate to real numbers and NUMVAL (normally ON) is ON. Variable treated as function of no arguments: E, PI. Functions of one argument: EXP, LOG, LN, LOG10, NORM, ARG, SQRT, RAD2DEG, RAD2DMS, DEG2RAD, DEG2DMS, DMS2DEG, DMS2RAD, SIN, ASIN, COS, ACOS, TAN, ATAN, COT, ACOT, SEC, ASEC, CSC, ACSC, SINH, ASINH, COSH, ACOSH, TANH, ATANH, COTH, ACOTH, SECH, ASECH, CSCH, ACSCH. Functions of two arguments: EXPT, LOGB, HYPOT, ATAN2. Function evaluation is carried out to the precision specified in the latest PRECISION statement. (The following functions are available only in symbolic mode, so they are not tested here: SIND, ASIND, COSD, ACOSD, TAND, ATAND, COTD, ACOTD, SECD, ASECD, CSCD, ACSCD, ATAN2D, CBRT); on rounded; precision 6; 12 a := exp 3; a := 20.0855 log a; 3.0 ln a; 3.0 log10 1000; 3.0 norm (-12.345); 12.345 % For real x, this is equivalent to ABS x. arg (-12.345); 3.14159 % For real x, this -> if x<0 then pi else 0.0. sqrt 3; 1.73205 ws**2; 3.0 deg2rad 30; 0.523599 rad2deg ws; 30.0 a := deg2dms 12.345; a := {12,20,42.0} % a will be a list. dms2deg ws; 12.345 dms2rad a; 0.215461 rad2deg ws; 12.345 asin 0.5; 0.523599 sin ws; 0.5 acos 0.5; 1.0472 cos ws; 0.5 atan 0.5; 0.463648 tan ws; 0.5 acot 0.5; 1.10715 cot ws; 0.5 asec 3; 1.23096 sec ws; 3.0 acsc 3; 0.339837 csc ws; 3.0 asinh 0.5; 0.481212 sinh ws; 0.5 acosh 2; 1.31696 cosh ws; 2.0 atanh 0.5; 0.549306 tanh ws; 0.5 acoth 2; 0.549306 coth ws; 2.0 sech 1; 0.648054 asech ws; 1 csch 1; 0.850918 acsch ws; 1 expt(2,1.234); 2.35218 logb(ws,2); 1.234 hypot(3,4); 5.0 a := -3*pi/4; a := - 2.35619 % Any -pi<a<=pi should work. atan2(sin a,cos a); - 2.35619 ws - a; 0 % Should be 0. precision 20; 6 % Functions will be computed to 20 places. sin 1.5; 0.99749498660405443094 asin ws; 1.5 precision 50; 20 % Functions computed to 50 places. sin 1.5; 0.99749498660405443094172337114148732270665142592212 asin ws; 1.5 precision 6; 50 comment If argument or value are complex, functions are not computed when dmode is ROUNDED; $ sin(1+i); sin(i + 1) % Complex argument. asin 2; asin(2) % Value would be complex. on complex; *** Domain mode rounded changed to complex-rounded % Now complex arguments and complex results will be handled. comment Complex functions of one argument: EXP, LOG, NORM, ARG, SQRT, SIN, ASIN, COS, ACOS, TAN, ATAN, COT, ACOT, SEC, ASEC, CSC, ACSC, SINH, ASINH, COSH. ACOSH, TANH, ATANH, COTH, ACOTH, SECH, ASECH, CSCH, ACSCH. (The following functions are available only in symbolic mode, so they are not tested here: SIND, ASIND, COSD, ACOSD, TAND, ATAND, COTD, ACOTD, SECD, ASECD, CSCD, ACSCD.) Complex function of two variables: EXPT, LOGB, ATAN2; e**(pi*i); - 1 + 1.22465e-16*i % Should be -1 (except for computational error.) log(1+i); 0.346574 + 0.785398*i exp ws; 1 + i norm(5*exp(2i)); 5.0 arg(5*exp(2i)); 2.0 sqrt(1+i); 1.09868 + 0.45509*i ws**2; 1 + i asin 2; 1.5708 - 1.31696*i sin ws; 2.0 - 1.06058e-16*i acos 2; 1.31696*i cos ws; 2.0 atan(1+i); 1.01722 + 0.402359*i tan ws; 1 + i acot(1+i); 0.553574 - 0.402359*i cot ws; 1 + i asec 0.1; 2.99322*i sec ws; 0.1 acsc 0.1; 1.5708 - 2.99322*i csc ws; 0.1 + 6.09254e-18*i sinh(1+i); 0.634964 + 1.29846*i asinh ws; 1 + i cosh(1+i); 0.83373 + 0.988898*i acosh ws; 1 + i atanh 2; 0.549306 + 1.5708*i tanh ws; 2.0 + 1.83697e-16*i acoth 0.3; 0.30952 + 1.5708*i coth ws; 0.3 - 5.57214e-17*i asech(1-i); 0.530638 + 1.11852*i sech ws; 1 - i acsch(1-i); 0.530638 + 0.452278*i csch ws; 1 - i expt(1+i,1-i); 2.80788 + 1.31787*i logb(ws,1+i); 1 - i a := 1+i; a := 1 + i % Any a such that - pi < repart a <= pi should work. atan2(sin a,cos a); 1 + i ws - a; 0 % Should be 0. clear a; % Further math package tests. %********************************************************************* %** %** The math package will compute the floating point values of ** %** the usual elementary functions, namely: ** %** sin asin sind asind sinh asinh ** %** cos acos cosd acosd cosh acosh ** %** tan atan tand atand tanh atanh ** %** cot acot cotd acotd coth acoth ** %** sec asec secd asecd sech asech ** %** csc acsc cscd acscd csch acsch ** %** atan2 atan2d ** %** exp ln sqrt ** %** expt log cbrt ** %** logb hypot ** %** log10 floor ** %** ceiling ** %** round ** %** ** %** All functions are computed to the accuracy of the floating- ** %** point precision of the system set up at the time. ** %** ** %********************************************************************* %** File #1===Trig Function Tests=== %** Trig functions are tested in both degrees and radians modes. %********************************************************************* symbolic; nil math!!label; "Math package mod 1.7, 1 May 93" symbolic procedure terpr(i,j); if remainder(i,j)=0 then terpri()$ on rounded; nil % We need !!plumin, etc. % #1: sind**2+cosd**2 test: ideal answers 1.0 for all args. for i:=0:45 do <<write " ",i,"->",sind float i**2+cosd float i**2; terpr(i,4)>>; 0->1.0 1->1.0 2->1.0 3->1.0 4->1.0 5->1.0 6->1.0 7->1.0 8->1.0 9->1.0 10->1.0 11->1.0 12->1.0 13->1.0 14->1.0 15->1.0 16->1.0 17->1.0 18->1.0 19->1.0 20->1.0 21->1.0 22->1.0 23->1.0 24->1.0 25->1.0 26->1.0 27->1.0 28->1.0 29->1.0 30->1.0 31->1.0 32->1.0 33->1.0 34->1.0 35->1.0 36->1.0 37->1.0 38->1.0 39->1.0 40->1.0 41->1.0 42->1.0 43->1.0 44->1.0 45->1.0 nil % #2: Quadrant test of sind, cosd: proper answers + +,+ -,- -,- +. begin scalar a; a:= sind 45.0; for i:= 0.0:3.0 do <<write " ",sind(i*90+45)/a," ", cosd (i*90+45)/a;terpri()>> end$ 1.0 1.0 1.0 -1.0 -1.0 -1.0 -1.0 1.0 % #3: Scaling test: all values should be 1 exactly. begin scalar a; a:= cosd 60.0; % for i:= -10.0:10.0 do write fix(cosd(60+i*360)/a)," " for i:= -10.0:10.0 do write round(cosd(60+i*360)/a)," " end$ 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 % #4: Test of radians -> degrees evaluation: ideal values 1.0. array a(6)$ begin for i:=1:6 do a(i):=sind(15.0*i); for i:=1:6 do <<write sin(!!pii2*i/6.0)/a(i)," "; terpr(i,3)>> end$ 1.0 1.0 1.0 1.0 1.0 1.0 % #5: Test of tand*cotd: ideal values 1.0. begin for i:=5 step 5 until 85 do <<write tand float i*cotd float i," "; terpr(i,25)>>; terpri() end$ 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 % #6: Test of secd**2-tand**2: ideal values 1.0. begin for i:=5 step 5 until 85 do <<write secd float i**2-tand float i**2," "; terpr(i,25)>> end$ 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 % #7: Test of cscd**2-cotd**2: ideal values 1.0. begin for i:=5 step 5 until 85 do <<write cscd float i**2-cotd float i**2," "; terpr(i,25)>> end$ 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 % #8: Test of asind+acosd: ideal values 1.0. begin write "sind and cosd"; terpri(); for i:=-10:10 do <<write (asind(0.1*i)+acosd(0.1*i))/90," "; terpr(i,5)>>; write "sin and cos";terpri(); for i:=-10:10 do <<write (acos(0.1*i)+asin(0.1*i))/!!pii2," "; terpr(i,5)>> end$ sind and cosd 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 sin and cos 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 % #9: Test of atand+acotd: ideal values 1.0. begin scalar x; write "tand, atand and acotd"; terpri(); for i:=-80 step 10 until 80 do <<x:=tand float i; write (atand x+acotd x)/90," "; terpr(i,50)>>; terpri(); write "tan, atan and acot";terpri(); for i:=-80 step 10 until 80 do <<x:= tan (!!pii2*i/90.0); write (atan x+acot x)/!!pii2," "; terpr(i,50)>> end$ tand, atand and acotd 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 tan, atan and acot 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 % #10: Test of atand tand: ideal values i for i:=-9:89. begin for i:=-9:89 do <<write " ",i,"->",if i=0 then 1.0 else atand tand float i; terpr(i,4)>> end$ -9->-9.0 -8->-8.0 -7->-7.0 -6->-6.0 -5->-5.0 -4->-4.0 -3->-3.0 -2->-2.0 -1->-1.0 0->1.0 1->1.0 2->2.0 3->3.0 4->4.0 5->5.0 6->6.0 7->7.0 8->8.0 9->9.0 10->10.0 11->11.0 12->12.0 13->13.0 14->14.0 15->15.0 16->16.0 17->17.0 18->18.0 19->19.0 20->20.0 21->21.0 22->22.0 23->23.0 24->24.0 25->25.0 26->26.0 27->27.0 28->28.0 29->29.0 30->30.0 31->31.0 32->32.0 33->33.0 34->34.0 35->35.0 36->36.0 37->37.0 38->38.0 39->39.0 40->40.0 41->41.0 42->42.0 43->43.0 44->44.0 45->45.0 46->46.0 47->47.0 48->48.0 49->49.0 50->50.0 51->51.0 52->52.0 53->53.0 54->54.0 55->55.0 56->56.0 57->57.0 58->58.0 59->59.0 60->60.0 61->61.0 62->62.0 63->63.0 64->64.0 65->65.0 66->66.0 67->67.0 68->68.0 69->69.0 70->70.0 71->71.0 72->72.0 73->73.0 74->74.0 75->75.0 76->76.0 77->77.0 78->78.0 79->79.0 80->80.0 81->81.0 82->82.0 83->83.0 84->84.0 85->85.0 86->86.0 87->87.0 88->88.0 89->89.0 % #11: Test of acot cotd: ideal values 10*i for i:=1:17. begin for i:=10 step 10 until 170 do <<write " ",i,"->",acotd cotd i; terpr(i,40)>>; terpri();terpri() end$ 10->10.0 20->20.0 30->30.0 40->40.0 50->50.0 60->60.0 70->70.0 80->80.0 90->90.0 100->100.0 110->110.0 120->120.0 130->130.0 140->140.0 150->150.0 160->160.0 170->170.0 % #12: Test of asind sind: ideal values 10*i for i:=-9:9. begin for i:=-90 step 10 until 90 do <<write " ",i,"->",asind sind float i; terpr(i,40)>> end$ -90->-90.0 -80->-80.0 -70->-70.0 -60->-60.0 -50->-50.0 -40->-40.0 -30->-30.0 -20->-20.0 -10->-10.0 0->0.0 10->10.0 20->20.0 30->30.0 40->40.0 50->50.0 60->60.0 70->70.0 80->80.0 90->90.0 % #13: Test of acosd cosd: ideal values 10*i for i:=1:18. begin for i:=10 step 10 until 180 do <<write " ",i,"->",acosd cosd float i; terpr(i,40)>> end$ 10->10.0 20->20.0 30->30.0 40->40.0 50->50.0 60->60.0 70->70.0 80->80.0 90->90.0 100->100.0 110->110.0 120->120.0 130->130.0 140->140.0 150->150.0 160->160.0 170->170.0 180->180.0 % #14: Test of acscd cscd: ideal values 10*i for i:=-9:9, except % error for i=0. begin for i:=-90 step 10 until 90 do <<write " ",i,"->",if i=0 then "error" else acscd cscd float i; terpr(i,40)>> end$ -90->-90.0 -80->-80.0 -70->-70.0 -60->-60.0 -50->-50.0 -40->-40.0 -30->-30.0 -20->-20.0 -10->-10.0 0->error 10->10.0 20->20.0 30->30.0 40->40.0 50->50.0 60->60.0 70->70.0 80->80.0 90->90.0 % #15: Test of asecd secd: ideal values 10*i for i :=0:18. except % error for i=9. begin for i:=0 step 10 until 180 do <<write" ",i,"->",if i=90 then "error" else asecd secd float i; terpr(i,40)>> end$ 0->0.0 10->10.0 20->20.0 30->30.0 40->40.0 50->50.0 60->60.0 70->70.0 80->80.0 90->error 100->100.0 110->110.0 120->120.0 130->130.0 140->140.0 150->150.0 160->160.0 170->170.0 180->180.0 %********************************************************************* %** ===Exp,Log,Sqrt,Cbrt, and Expt Function tests=== %********************************************************************* % #16: Test of properties of exp function: ideal results 1.0. array b(5)$ begin scalar x; x:=0; write "multiplicative property";terpri(); for i:=0:5 do b(i):=1+i/6.0; for i:=0:5 do for j:=i:5 do <<write " ",exp (b(i)+b(j))/exp(b(i))/exp(b(j)); terpr(x:=x+1,5)>> end$ multiplicative property 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 % #17: Various properties of exp: ideal results 1.0. begin write "inverse property"$ terpri()$ for i:=1:5 do write " ",exp(b(i))*exp(-b(i));terpri(); write "squares"; terpri(); for i:=-10:10 do <<write " ",sqrt(exp(0.2*i))/exp(0.1*i); terpr(i,5)>>; write "cubes"; terpri(); for i:=-10:10 do <<write " ",cbrt(exp(0.3*i))/exp(0.1*i); terpr(i,5)>> end$ inverse property 1.0 1.0 1.0 1.0 1.0 squares 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 cubes 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 % #18: Test of log exp: ideal results 1.0. begin for i:=-5:5 do <<write if i=0 then "0/0" else (log exp float i)/i," "; terpr(i,5)>> end$ 1.0 1.0 1.0 1.0 1.0 0/0 1.0 1.0 1.0 1.0 1.0 % #19: Test of log10 expt(10.0,i): ideal results 1.0. begin scalar i; write "small values i:=-5:5"; terpri(); for j:=-5:5 do <<write if j neq 0 then log10 float expt(10.0,j)/j else "zero"," "; terpr(j,5)>>; write "large i=2**j where j:=0:6"; terpri(); for j:=0:5 do <<write (log10 float expt(10.0,2**j))/2**j," "; terpr(j,5)>>; terpri(); write "noninteger values of i=j/10.0 where j:=1:20";terpri(); for j:=1:20 do <<i:=j/10.0; write (log10 float expt(10,i))/i," "; terpr(j,5)>> end$ small values i:=-5:5 1.0 1.0 1.0 1.0 1.0 zero 1.0 1.0 1.0 1.0 1.0 large i=2**j where j:=0:6 1.0 1.0 1.0 1.0 1.0 1.0 noninteger values of i=j/10.0 where j:=1:20 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 % #20: Test of properties of expt(x,i)*(expt(x,-i). ideal result 1.0. begin integer j; for x:=2:6 do for i:=2:6 do <<write expt(float x,i)*expt(float x,-i)," "; terpr(j:=j+1,5)>> end$ 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 % #21: Test of expt(-x,i)/expt(x,i) for fractional i. begin integer j,k; write "odd numerator. ideal result -1.0"; terpri(); for i:=1:10 do <<k:=(2*i-1.0)/(8*i+1); write rexpt(-8,k)/rexpt(8,k)," "; terpr(j:=j+1,5)>>; write "even numerator. ideal result 1.0"; terpri(); for i:=1:10 do <<k:=(2.0*i)/(8*i+1); write rexpt(-8,k)/rexpt(8,k)," "; terpr(j:=j+1,5)>> end$ odd numerator. ideal result -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 even numerator. ideal result 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 % #22: Test of properties of ln or log or logb: % inverse argument: ideal result -1.0. begin integer x; for i:=2:5 do for j:= 2:10 do <<x:=x+1; write logb(float i,float j)/logb(float i,1.0/j)," "; terpr(x,5)>> end$ -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 -1.0 % #23: Test of log(a*b) = log a+log b: ideal result 1.0. begin integer x; for i:=1:5 do for j:=i:5 do <<write log (i*j*0.01)/(log (i*0.1)+log(j*0.1))," "; terpr(x:=x+1,5)>> end$ 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 % #24: Test of sqrt x*sqrt x/x for x:=5i*(5i/3)**i where i:=1:20 % (test values strictly arbitrary): ideal results 1.0. begin scalar x,s; for i:=1:20 do <<x:= 5.0*i;s:=sqrt(x:=x*(expt(x/3,i))); write s*s/x," "; terpr(i,5)>> end$ 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 % #25: Test of cbrt x**3/x for x:=5i*(5i/3)**i where i:=-9:10 % (test values strictly arbitrary):ideal results 1.0. begin scalar x,s; for i:=-9:10 do <<x:= 5.0*i; if i neq 0 then s:= cbrt(x:=x*(expt(x/3,i))); write if i=0 then 1 else s*s*s/x," "; terpr(i,5)>> end$ 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 %********************************************************************* %** ===Hyperbolic Function Tests=== %********************************************************************* % #26: Test of sinh x+ cosh x= exp x: ideal results 1.0. begin scalar x; for i:=1:10 do <<x:=ln float i$ write (sinh x+cosh x)/exp x," "$ terpr(i,5)>> end$ 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 % #27: Test of cosh x-sinh x= exp(-x): ideal results 1.0. begin scalar x; for i:=1:10 do <<x:=ln float i$ write(cosh x-sinh x)*exp x," "$ terpr(i,5)>> end$ 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 % #28: Test of (cosh x)**2-(sinh x)**2: ideal results 1.0. begin scalar x$ for i:=1:10 do <<x:=ln float i$write(cosh x)**2-(sinh x)**2," "; terpr(i,5)>> end$ 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 % #29: Test of tanh*cosh/sinh: ideal results 1.0. begin scalar x; for i:=1:20 do <<x:=ln(i*0.1); write if i=10 then 1 else tanh x*cosh x/sinh x," "; terpr(i,5)>> end$ 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 % #30: Test of tanh*coth: ideal results 1.0. begin scalar x; for i:=1:20 do <<x:=ln(i*0.1); write if i=10 then 1 else tanh x*coth x," "; terpr(i,5)>> end$ 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 % #31: Test of sech*cosh: ideal results 1.0. begin scalar x; for i:=1:20 do <<x:=ln(i*0.1); write sech x*cosh x," "; terpr(i,5)>> end$ 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 % #32: Test of csch*sinh: ideal results 1.0. begin scalar x; for i:=1:20 do <<x:=ln(i*0.1); write if i=10 then 1 else csch x*sinh x," "; terpr(i,5)>> end$ 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 % #33: Test of asinh sinh: ideal results 1.0. begin scalar x; for i:=1:20 do <<x:=ln(i*0.1); write if i=10 then 1 else asinh sinh x/x," "; terpr(i,5)>> end$ 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 % #34: Test of acosh cosh: ideal results 1.0. However, acosh x % loses accuracy as x -> 1 since d/dx cosh x -> 0. begin scalar x; for i:=1:20 do <<x:=ln(1+i*0.05); write acosh cosh x/x," "; terpr(i,5)>> end$ 1.0 1.0 0.99999999999999 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 % #35: Test of cosh acosh:ideal results 1.0. begin scalar x; for i:=1:50 do <<x:=1+i/25.0; write (cosh acosh x)/x," "; terpr(i,5)>> end$ 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 % #36: Test of atanh tanh: ideal results 1.0. begin scalar x; for i:=1:20 do <<x:=ln(i*0.1); write if i=10 then 1 else (atanh tanh x)/x," "; terpr(i,5)>> end$ 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 % #37: Test of acoth coth: ideal results 1.0. begin scalar x; for i:=1:20 do <<x:=ln(i*0.1); write if i=10 then 1 else (acoth coth x)/x," "; terpr(i,5)>> end$ 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 % #38: Test of asech sech: ideal results 1.0. However, asech x % loses accuracy as x -> 1 since d/dx sech x -> 0. begin scalar x; for i:=1:20 do <<x:=ln(1+i*0.05); write (asech sech x)/x," "; terpr(i,5)>> end$ 1.0 1.0 0.99999999999999 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 % #39: Test of acsch csch: ideal results 1.0. begin scalar x; for i:=1:20 do <<x:=ln(i*0.1); write if i=10 then 1 else (acsch csch x)/x," "; terpr(i,5)>> end$ 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 end; Time for test: 50 ms, plus GC time: 10 ms