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
