From 98cab31fc3659e33aef260efca55bf9f1753164c Mon Sep 17 00:00:00 2001
From: Lars-Dominik Braun <lars@6xq.net>
Date: Mon, 11 Feb 2019 11:49:19 +0100
Subject: Add source files from Michael

---
 system/base/unknown/src/mathlib | 359 ++++++++++++++++++++++++++++++++++++++++
 1 file changed, 359 insertions(+)
 create mode 100644 system/base/unknown/src/mathlib

(limited to 'system/base/unknown/src/mathlib')

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; 
-- 
cgit v1.2.3