(* ------------------- VERSION 2     06.03.86 ------------------- *)
PACKET mathlib DEFINES sqrt, **, exp, ln, log2, log10, e, pi,
                       sin, cos, tan, sind, cosd, tand,
                       arctan, arctand, random, initializerandom :
 
LET  pii   =  3.141592653589793238462,
     pi2   =  1.570796326794896619231,
     pi3   =  1.047197551196597746154,
     pi6   =  0.523598775598298873077,
     pi4   =  1.273239544735162686151,
     ln2   =  0.693147180559945309417,
     lg2   =  0.301029995663981195213,
     ln10  =  2.302585092994045684018,
     lge   =  0.434294481903251827651,
     ei    =  2.718281828459045235360,
     pi180 = 57.295779513082320876798,
     sqrt3 =  1.732050807568877293527,
     sqr3  =  0.577350269189625764509,
     sqr3p2=  3.732050807568877293527,
     sqr3m2=  0.267949192431122706473,
     sqr2  =  0.707106781186547524400;

REAL VAR rdg::0.4711; 
 
REAL PROC pi:   pii     END PROC pi; 
REAL PROC e :   ei      END PROC e; 
 
REAL PROC  ln  ( REAL CONST x ): 
     log2(x) * ln2
END PROC ln; 
 
REAL PROC log10( REAL CONST x ): 
     log2(x) * lg2 
END PROC log10; 

REAL PROC log2 ( REAL CONST z ): 
 REAL VAR t, summe::0.0, x::z;
 IF   x=1.0  THEN 0.0
 ELIF x>0.0  THEN normal
             ELSE errorstop("log2: " + text (x,20)); 0.0   FI. 

normal: 
  IF   x >= 0.5   THEN  normalise downwards 
                  ELSE  normalise upwards      FI; 
  IF   x < sqr2   THEN  summe := summe - 0.75; t := trans8
                  ELSE  summe := summe - 0.25; t := trans2  FI;
  summe + reihenentwicklung. 
 
 normalise downwards: 
  WHILE  x >= 8.0 REP    x := 0.0625 * x;  summe:=summe+4.0  PER; 
  WHILE  x >= 1.0 REP    x :=    0.5 * x;  summe:=summe+1.0  PER. 
 
 normalise upwards: 
  WHILE x<=0.0625 REP    x :=   16.0 * x;  summe:=summe-4.0  PER; 
  WHILE x<= 0.5   REP    x :=    2.0 * x;  summe:=summe-1.0  PER. 
 
 trans8: (x - 0.5946035575013605)/(x + 0.5946035575013605). 
 trans2: (x - 0.8408964152537145)/(x + 0.8408964152537145).
 
 reihenentwicklung:  x := t * t;         t * 0.06405572387119384648 * 
  ((((((3.465*x+4.095)*x+5.005)*x+6.435)*x+9.009)*x+15.015)*x+45.045)
END PROC log2;
 
REAL PROC sqrt ( REAL CONST z ): 
  REAL VAR y0, y1, x::z; 
  INT  VAR p :: decimal exponent(x) DIV 2; 
  IF   p <= -64  THEN 0.0
  ELIF x < 0.0  THEN errorstop("sqrt: " + text (x,20)); 0.0
                ELSE nontrivial    FI. 
 
 nontrivial: 
  set exp (decimal exponent (x) -p-p, x);
  IF x<10.0 THEN x := 5.3176703 - 40.760905/( 8.408065 + x )
            ELSE x := 16.81595  - 1288.973 /( 84.08065 + x )  FI; 
  y0 := x;
  set exp (decimal exponent (x) + p, y0);
  y1 := 0.5 * ( y0 + z/y0 ); 
  y0 := 0.5 * ( y1 + z/y1 ); 
  y1 := 0.5 * ( y0 + z/y0 ); 
        0.5 * ( y1 + z/y1 ) 
END PROC sqrt; 
 
REAL PROC exp ( REAL CONST z ): 
 REAL VAR x::z, a::1.0;  BOOL VAR negativ :: x<0.0;
 IF   negativ        THEN x := -x  FI;
 IF   x>292.42830676 
   THEN IF NOT negativ THEN errorstop ("REAL-Ueberlauf") FI ; 0.0
 ELIF x<=0.0001 
   THEN ( 0.5*z + 1.0 ) * z + 1.0
   ELSE approx 
 FI.
 
 approx: 
 IF x > ln10 
 THEN x := lge*x; 
      a := 1.0;
      set exp (int(x), a);
      x := frac(x)*ln10
 FI;
 IF x >= 2.0   THEN a := 7.389056098930650227230*a; x := x-2.0   FI;
 IF x >= 1.0   THEN a := 2.718281828459045235360*a; x := x-1.0   FI;
 IF x >= 0.5   THEN a := 1.648721270700128146848*a; x := x-0.5   FI;
 IF x >= 0.25  THEN a := 1.284025416687741484073*a; x := x-0.25  FI;
 IF x >= 0.125 THEN a := 1.133148453066826316829*a; x := x-0.125 FI;
 IF x >= 0.0625THEN a := 1.064494458917859429563*a; x := x-0.0625FI;
 a:=a/50.4*(((((((0.01*x+0.07)*x+0.42)*x+2.1)*x+8.4)*x+25.2)*x+50.4)*x+50.4);
 IF negativ THEN 1.0/a ELSE a FI .

ENDPROC exp ;

REAL PROC tan (REAL CONST x): 
  IF x < 0.0 THEN - tg( -x * pi4)
             ELSE   tg(  x * pi4)  FI
END PROC tan; 
 
REAL PROC tand (REAL CONST x): 
  IF x < 0.0 THEN - tg( -x / 45.0)
             ELSE   tg(  x / 45.0) FI
END PROC tand; 
 
REAL PROC tg (REAL CONST x ): 
  REAL VAR q::floor(x), s::x-q; INT VAR n;
  q := q - floor(0.25*q) * 4.0 ; 
  IF q < 2.0
  THEN IF q < 1.0
       THEN n:=0;
       ELSE n:=1; s := 1.0 - s FI
  ELSE IF q < 3.0
       THEN n:=2;
       ELSE n:=3; s := 1.0 - s FI
  FI;
  q := s * s;
  q := (((((((((-5.116186989653120e-11*q-5.608325022830701e-10)*q- 
  9.526170109403018e-9)*q-1.517906721393745e-7)*q-2.430939946375515e-6)*q- 
  3.901461426385464e-5)*q-6.324811612385572e-4)*q-1.076606829172646e-2)*q- 
  0.2617993877991508)*q+pi4); 

  SELECT n OF
  CASE 0 : s/q
  CASE 1 : q/s
  CASE 2 : -q/s
  OTHERWISE : -s/q ENDSELECT  .

END PROC tg;

REAL PROC sin ( REAL CONST x ): 
  REAL VAR y, r, q;
  IF x < 0.0 THEN y := -x; q := 4.0 ELSE y := x; q := 0.0 FI;
  y := y * pi4;
  r := floor(y);
  sincos( q+r , y-r )
END PROC sin; 
 
REAL PROC sind ( REAL CONST x ): 
  REAL VAR y, r, q; 
  IF x < 0.0 THEN y := -x; q := 4.0 ELSE y := x; q := 0.0 FI;
  y := y / 45.0;
  r := floor(y);
  sincos( q+r , y-r )
END PROC sind; 
 
REAL PROC cos ( REAL CONST x ): 
  REAL VAR y, q; 
  IF x < 0.0 THEN y := -x ELSE y := x FI;
  y := y * pi4;
  q := floor(y); 
  sincos( q+2.0, y-q )
END PROC cos; 
 
REAL PROC cosd ( REAL CONST x ): 
  REAL VAR y, q; 
  IF x < 0.0 THEN y := -x ELSE y := x FI;
  y := y / 45.0;
  q := floor(y);
  sincos( q+2.0, y-q )
END PROC cosd; 
 
REAL PROC sincos ( REAL CONST q, y ): 
 REAL VAR r :: q - floor( 0.125*q + 0.1 ) * 8.0;
 IF r >= 4.0 THEN IF r >= 6.0 THEN IF r >= 7.0 THEN - sin approx(1.0-y)
                                               ELSE - cos approx(y)  FI
                              ELSE IF r >= 5.0 THEN - cos approx(1.0-y)
                                               ELSE - sin approx(y)  FI FI
             ELSE IF r >= 2.0 THEN IF r >= 3.0 THEN   sin approx(1.0-y)
                                               ELSE   cos approx(y)  FI
                              ELSE IF r >= 1.0 THEN   cos approx(1.0-y)
                                               ELSE   sin approx(y)  FI FI FI
END PROC sincos; 
 
REAL PROC sin approx ( REAL CONST x ): 
 REAL VAR z::x*x; 
 x*((((((0.6877101540593035e-11*z-0.1757149296873372e-8)*z+0.3133616216672568
 e-6)*z-0.3657620415845891e-4)*z+0.2490394570188737e-2)*z-0.807455121882e-1)*
 z+0.7853981633974483)
END PROC sin approx; 
 
REAL PROC cos approx ( REAL CONST x ): 
 REAL VAR z::x*x; 
 ((((((-0.3857761864560276e-12*z+0.115004970178141e-9)*z-0.246113638267419e-7
 )*z+0.3590860445885748e-5)*z-0.3259918869266875e-3)*z+0.1585434424381541e-1)
 *z-0.3084251375340425)*z+1.0
END PROC cos approx; 
 
REAL PROC arctan ( REAL CONST y ): 
 REAL VAR f, z, x; BOOL VAR neg :: y < 0.0;
 IF   neg   THEN x := -y ELSE x := y FI;
 IF   x>1.0 THEN f := a ELSE f := -b; neg := NOT neg FI;
 z := x * x;
 x := x/(((((((0.0107090276046822*z-0.01647757182108040)*z 
      +0.02177846332482151)*z-0.03019339673273880)*z+0.04656083561183398)*z 
      -0.0888888888888888)*z+0.3333333333333333)*z+1.0);
 IF   neg   THEN  x - f  ELSE  f - x  FI.

 a:IF x>sqr3p2 THEN x := 1.0/x; pi2 ELSE x := 4.0/(sqrt3+x+x+x)-sqr3; pi3 FI.
 b:IF x<sqr3m2 THEN   0.0           ELSE x := sqrt3 - 4.0/(sqrt3+x);  pi6 FI
END PROC arctan; 
 
REAL PROC arctand ( REAL CONST x ): 
  arctan(x) * pi180
END PROC arctand; 
 
REAL OP ** ( REAL CONST b, e ): 
 IF b=0.0 
   THEN IF e=0.0 THEN 1.0 ELSE 0.0 FI
 ELIF b < 0.0 
   THEN errorstop("("+text(b,20)+") ** "+text(e)); (-b) ** e
 ELSE   exp( e * log2( b ) * ln2 )
 FI
END OP **;
 
REAL OP ** ( REAL CONST a, INT CONST b ) : 

  REAL VAR p := 1.0 ,
           r := a ;
  INT VAR  n := ABS b ,
           m ; 
  IF (a = 0.0 OR a = -0.0) 
    THEN IF b = 0
           THEN 1.0
           ELSE 0.0
         FI
    ELSE WHILE n>0 REP
           m := n DIV 2 ;
           IF m + m = n
             THEN n := m ;
                  r := r*r 
             ELSE n DECR 1 ;
                  p := p*r
           FI
         END REP ; 
         IF b>0
           THEN p
           ELSE 1.0 / p
         FI
  FI .

END OP ** ; 
 
REAL PROC random: 
  rdg:=rdg+pii;rdg:=rdg*rdg;rdg:=rdg*rdg;rdg:=rdg*rdg;rdg:=frac(rdg);rdg
END PROC random; 
 
PROC initializerandom ( REAL CONST z ): 
  rdg := frac(z) 
END PROC initializerandom; 
 
END PACKET mathlib;