Initial import

1 files changed, 268 insertions, 0 deletions

diff --git a/system/base/1.7.5/src/mathlib b/system/base/1.7.5/src/mathlib
new file mode 100644
index 0000000..c726495
--- /dev/null
+++ b/system/base/1.7.5/src/mathlib
@@ -0,0 +1,268 @@
+(* ------------------- 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; 
+