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;