Add source files from Michael

1 files changed, 359 insertions, 0 deletions

diff --git a/system/base/unknown/src/mathlib b/system/base/unknown/src/mathlib
new file mode 100644
index 0000000..be44ff6
--- /dev/null
+++ b/system/base/unknown/src/mathlib
@@ -0,0 +1,359 @@
+ 

+PACKET mathlib DEFINES sqrt,**,exp,ln,log2,log10,sin,cos, 

+                       tan,arctan,sind,cosd,tand,arctand,e,pi,

+                       random,initializerandom :

+ 

+ 

+REAL VAR rdg::0.4711; 

+ 

+REAL PROC pi: 

+   3.141592653589793. 

+END PROC pi; 

+ 

+REAL PROC e: 

+  2.718281828459045. 

+END PROC e; 

+ 

+REAL PROC ln(REAL CONST x): 

+LET ln2= 0.6931471805599453; 

+log2(x)*ln2. 

+END PROC ln; 

+ 

+REAL PROC log2(REAL CONST z): 

+INT VAR k::0,p::0; 

+REAL VAR m::0.0,x::z,t::0.0,summe::0.0; 

+IF   x>0.0 

+THEN normal 

+ELSE errorstop("log2 mit negativer zahl");4711.4711 

+FI. 

+normal: 

+ IF   x>=0.5 

+ THEN normalise downwards 

+ ELSE normalise upwards 

+  FI; 

+  IF   x>=0.1 AND x< 0.7071067811865475 THEN 

+       t:=(x-0.5946035575013605)/(x+0.5946035575013605); 

+       summe:=reihenentwicklung (t) - 0.75 

+  FI; 

+  IF  x>=0.7071067811865475 AND x < 1.0 THEN 

+      t:=(x - 0.8408964152537145)/(x+0.8408964152537145); 

+      summe:= reihenentwicklung(t)-0.25 

+  FI; 

+  summe-real(p - 4*k). 

+ 

+ normalise downwards: 

+  WHILE  x>= 16.0 REP 

+    x:=x/16.0;k:=k+1; 

+  END REP; 

+  WHILE  x>=0.5 REP 

+    x:=x/2.0;p:=p-1; 

+  END REP. 

+ 

+ normalise upwards: 

+  WHILE x<=0.0625 REP 

+    x:=x*16.0;k:=k-1; 

+  END REP; 

+  WHILE x<= 0.5 REP 

+    x:=x*2.0;p:=p+1; 

+  END REP. 

+ 

+END PROC log2; 

+ 

+REAL PROC reihenentwicklung(REAL CONST x):

+  REAL VAR i::39.0,s::1.0/39.0; 

+  LET ln2=0.6931471805599453; 

+  WHILE i>1.0 REP 

+    i:=i-2.0;s:=s*x*x + 1.0/i; 

+  END REP; 

+  s*2.0*x/ln2. 

+END PROC reihenentwicklung; 

+ 

+REAL PROC log10(REAL CONST x): 

+  LET lg2=0.301029995664; 

+  log2(x)*lg2. 

+END PROC log10; 

+ 

+REAL PROC sqrt(REAL CONST z): 

+  REAL VAR y0,y1,x::z; 

+  INT VAR p::0; 

+  BOOL VAR q::FALSE; 

+  IF   x<0.0 

+  THEN errorstop("sqrt von negativer zahl");0.0

+  ELSE correct 

+  FI. 

+ 

+ correct: 

+  IF   x=0.0 

+  THEN 0.0 

+  ELSE nontrivial 

+  FI. 

+ 

+ nontrivial: 

+  IF   x<0.01 

+  THEN small 

+  ELSE notsmall 

+  FI. 

+ 

+ 

+ notsmall: 

+  IF   x>1.0 

+  THEN big 

+  ELSE normal 

+  FI. 

+ 

+ small: 

+  WHILE x<0.01 REP 

+    p:=p-1;x:=x*100.0; 

+  END REP; 

+  normal. 

+ 

+ big: 

+  WHILE x>=1.0 REP

+    p:=p+1;x:=x/100.0; 

+  END REP; 

+  normal. 

+ 

+ normal: 

+  IF   x<0.1 

+  THEN x:=x*10.0;q:=TRUE 

+  FI; 

+  y0:=10.0**p*(1.681595-1.288973/(0.8408065+x)); 

+  IF   q 

+  THEN y0:=y0/3.162278 

+  FI; 

+  y1:=(y0+z/y0)/2.0; 

+  y0:=(y1+z/y1)/2.0; 

+  y1:=(y0+z/y0)/2.0; 

+  (y1-z/y1)/2.0+z/y1. 

+ 

+END PROC sqrt; 

+ 

+REAL PROC exp(REAL CONST z): 

+  REAL VAR c,d,x::z, a, b ; 

+  IF   x<-180.2187 

+  THEN 0.0 

+  ELIF x<0.0 

+  THEN 1.0/exp(-x) 

+  ELIF x=0.0 

+  THEN 1.0 

+  ELSE x:=x/0.6931471805599453;approx 

+  FI. 

+ 

+ approx: 

+  a:=floor(x/4.0)+1.0; 

+  b:=floor(4.0*a-x); 

+  c:=(4.0*a-b-x)*16.0; 

+  d:=(c -floor(c))/16.0; 

+  d:=d*0.6931471805599453; 

+  ( (16.0 POWER a) / (2.0 POWER b) / (1.044273782427419 POWER c ))* 

+  ((((((0.135910788320380e-2*d-0.8331563191293753e-2)*d 

+  +0.4166661437490328e-1)*d-0.1666666658727157)*d+0.4999999999942539)*d 

+  - 0.9999999999999844)*d+1.0). 

+ 

+ENDPROC exp ; 

+ 

+REAL OP POWER (REAL CONST basis, exponent) : 

+ 

+  IF floor (exponent) = 0.0 

+    THEN 1.0 

+    ELSE power 

+  FI . 

+ 

+power : 

+  REAL VAR counter := floor (abs (exponent)) - 1.0 , 

+           result  := basis ; 

+  WHILE counter > 0.0 REP 

+    result := result * basis ; 

+    counter := counter - 1.0 

+  PER ; 

+  IF exponent > 0.0 

+    THEN result 

+    ELSE 1.0 / result 

+  FI . 

+ 

+ENDOP POWER ; 

+ 

+REAL PROC tan (REAL CONST x): 

+  REAL VAR p; 

+  p:=1.273239544735168*ABSx; 

+  tg(p)*sign(x). 

+END PROC tan; 

+ 

+REAL PROC tand(REAL CONST x): 

+  REAL VAR p; 

+  p:=0.02222222222222222*ABSx; 

+  tg(p)*sign(x). 

+END PROC tand; 

+ 

+REAL PROC tg(REAL CONST x): 

+  REAL VAR r,s,u,q; 

+  q:=floor(x);r:=x-q; 

+  IF q = floor(q/2.0) * 2.0 

+  THEN s:=r 

+  ELSE s:=(1.0-r) 

+  FI; 

+  q:= q - floor(q/4.0) * 4.0 ; 

+  u:=s*s; 

+  s:=s*0.785398163397448; 

+  s:=s/(((((((((-0.4018243865271481e-10*u-0.4404768172667185e-9)*u- 

+  0.748183650813680e-8)*u-0.119216115119129e-6)*u-0.1909255769212821e-5)*u- 

+0.3064200638849133e-4)*u-0.4967495424202482e-3)*u-0.8455650263333471e-2)*u- 

+  0.2056167583560294)*u+1.0); 

+  IF   q=0.0 

+  THEN s 

+  ELIF q=3.0 

+  THEN -s 

+  ELIF q=1.0 

+  THEN 1.0/s 

+  ELSE -1.0/s 

+  FI . 

+ 

+END PROC tg;

+ 

+REAL PROC sin(REAL CONST x): 

+  REAL VAR y,r; 

+  INT VAR q; 

+  y:=ABS x*1.273239544735168; 

+  q:=int(y); 

+  r:=y-real(q); 

+  IF   x<0.0 

+  THEN q:=q+4 

+  FI; 

+  sincos(q,r). 

+END PROC sin; 

+ 

+REAL PROC sind(REAL CONST x): 

+  REAL VAR y,r; 

+  INT VAR q; 

+  y:=ABSx/45.0; 

+  q:=int(y); 

+  r:=y-real(q); 

+  IF   x<0.0 

+  THEN q:=q+4 

+  FI; 

+  sincos(q,r). 

+END PROC sind; 

+ 

+ 

+REAL PROC cos(REAL CONST x): 

+  REAL VAR y,r; 

+  INT VAR q; 

+  y:=ABS x*1.273239544735168; 

+  q:=int(y); 

+  r:=y-real(q); 

+  q:=q+2; 

+  sincos(q,r). 

+END PROC cos; 

+ 

+REAL PROC cosd(REAL CONST x): 

+  REAL VAR y,r; 

+  INT VAR q; 

+  y:=ABS x/45.0; 

+  q:=int(y); 

+  r:=y-real(q); 

+  q:=q+2; 

+  sincos(q,r). 

+END PROC cosd; 

+ 

+ 

+REAL PROC sincos(INT VAR q,REAL VAR r): 

+  SELECT q MOD 8 + 1 OF 

+  CASE 1 : sin approx(r) 

+  CASE 2 : cos approx (1.0-r) 

+  CASE 3 : cos approx(r) 

+  CASE 4 : sin approx(1.0-r) 

+  CASE 5 : - sin approx(r) 

+  CASE 6 : - cos approx(1.0-r) 

+  CASE 7 : - cos approx(r) 

+  CASE 8 : - sin approx(1.0-r) 

+  OTHERWISE 0.0 

+  END SELECT 

+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.313361621667256

+8

+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 x): 

+REAL VAR z::x*x; 

+IF   x<0.0 THEN -arctan(-x) 

+ELIF x>1.0 THEN 3.141592653589792/2.0-arctan(1.0/x) 

+ELIF x*1.0e16>2.67949192431e15   THEN pi/6.0+arctan(1.732050807568877-4.0 

+/(x+1.732050807568877)) 

+ELSE  x/(((((((0.0107090276046822*z-0.01647757182108040)*z 

+      +0.02177846332482151)*z-0.03019339673273880)*z+0.04656083561183398)*z 

+      -0.0888888888888888)*z+0.3333333333333333)*z+1.0)FI. 

+END PROC arctan; 

+ 

+REAL PROC arctand(REAL CONST x): 

+  arctan(x)/3.1415926589793*180.0. 

+END PROC arctand; 

+ 

+ 

+BOOL PROC even(INT CONST number): 

+ (number DIV 2)*2=number. 

+END PROC even; 

+ 

+REAL OP **(REAL CONST base,exponent): 

+  IF   base<0.0 

+  THEN errorstop("hoch mit negativer basis") 

+  FI; 

+  IF   base=0.0 

+  THEN test exponent 

+  ELSE

+      exp(exponent*ln(base)) 

+  FI. 

+ 

+ test exponent: 

+  IF   exponent=0.0 

+  THEN errorstop("0**0 geht nicht");4711.4711 

+  ELSE 0.0 

+  FI. 

+ 

+END OP **; 

+ 

+ 

+REAL PROC sign(REAL CONST number): 

+  IF   number >0.0  THEN 1.0 

+  ELIF number <0.0  THEN -1.0 

+  ELSE 0.0 

+  FI. 

+END PROC sign ;

+ 

+REAL OP **(REAL CONST a,INT CONST b): 

+REAL VAR p::1.0,r::a;INT VAR n::ABS b; 

+WHILE n>0 REP 

+  IF   n MOD 2=0 

+  THEN n:=n DIV 2;r:=r*r 

+  ELSE n DECR 1;p:=p*r 

+  FI; 

+END REP; 

+IF   b>0 

+THEN p 

+ELSE 1.0/p 

+FI. 

+END OP **; 

+ 

+ 

+ 

+REAL PROC random: 

+rdg:=rdg+pi;rdg:=rdg*rdg;rdg:=rdg*rdg;rdg:=rdg*rdg;rdg:=frac(rdg);rdg. 

+END PROC random; 

+ 

+ 

+PROC initializerandom(REAL CONST z): 

+ rdg:=z; 

+END PROC initializerandom; 

+ 

+END PACKET mathlib;