d0p1/ack
1
0
Fork 0
Code Issues Pull requests Projects Releases Packages Wiki Activity

Used new math lib of C to create new version of Mathlib

Browse source
This commit is contained in:
ceriel 1988-07-25 16:41:51 +00:00
parent e98a670850
commit dbbff76f4c
2 changed files with 253 additions and 127 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

30
lang/m2/libm2/Mathlib.def
View file

@ -8,19 +8,27 @@ DEFINITION MODULE Mathlib;
(* Some mathematical constants: *) (* Some mathematical constants: *)
CONST CONST
(* From: Handbook of Mathematical Functions (* From: Computer Approximations
Edited by M. Abramowitz and I.A. Stegun Hart, Cheney, e.a.
National Bureau of Standards The SIAM Series in Applied Mathematics
Applied Mathematics Series 55 John Wiley & Sons, INC. New York London Sydney, 1968
*) *)
pi = 3.141592653589793238462643; pi = 3.14159265358979323846264338327950288;
twicepi = 6.283185307179586476925286; twicepi = 6.28318530717958647692528676655900576;
halfpi = 1.570796326794896619231322; halfpi = 1.57079632679489661923132169163975144;
quartpi = 0.785398163397448309615661; quartpi = 0.78539816339744830961566084581987572;
e = 2.718281828459045235360287; e = 2.71828182845904523536028747135266250;
ln2 = 0.693147180559945309417232; ln2 = 0.69314718055994530941723212145817657;
ln10 = 2.302585092994045684017992; ln10 = 2.30258509299404568401799145468436421;
longpi = 3.14159265358979323846264338327950288D;
longtwicepi = 6.28318530717958647692528676655900576D;
longhalfpi = 1.57079632679489661923132169163975144D;
longquartpi = 0.78539816339744830961566084581987572D;
longe = 2.71828182845904523536028747135266250D;
longln2 = 0.69314718055994530941723212145817657D;
longln10 = 2.30258509299404568401799145468436421D;
(* basic functions *) (* basic functions *)

350
lang/m2/libm2/Mathlib.mod
View file

@ -14,15 +14,11 @@ IMPLEMENTATION MODULE Mathlib;
FROM EM IMPORT FIF, FEF; FROM EM IMPORT FIF, FEF;
FROM Traps IMPORT Message; FROM Traps IMPORT Message;
(* From: Handbook of Mathematical Functions
Edited by M. Abramowitz and I.A. Stegun
National Bureau of Standards
Applied Mathematics Series 55
*)
CONST CONST
OneRadianInDegrees = 57.295779513082320876798155D; OneRadianInDegrees = 57.295779513082320876798155D;
OneDegreeInRadians = 0.017453292519943295769237D; OneDegreeInRadians = 0.017453292519943295769237D;
Sqrt2 = 1.41421356237309504880168872420969808D;
OneOverSqrt2 = 0.70710678118654752440084436210484904D;
(* basic functions *) (* basic functions *)
@ -32,8 +28,7 @@ IMPLEMENTATION MODULE Mathlib;
END pow; END pow;
PROCEDURE longpow(x: LONGREAL; i: INTEGER): LONGREAL; PROCEDURE longpow(x: LONGREAL; i: INTEGER): LONGREAL;
VAR VAR val: LONGREAL;
val: LONGREAL;
ri: LONGREAL; ri: LONGREAL;
BEGIN BEGIN
ri := FLOATD(i); ri := FLOATD(i);
@ -93,7 +88,7 @@ IMPLEMENTATION MODULE Mathlib;
temp := temp / 2.0D; temp := temp / 2.0D;
exp := exp + 2; exp := exp + 2;
END; END;
FOR i := 0 TO 4 DO FOR i := 0 TO 5 DO
temp := 0.5D*(temp + x/temp); temp := 0.5D*(temp + x/temp);
END; END;
RETURN temp; RETURN temp;
@ -105,42 +100,50 @@ IMPLEMENTATION MODULE Mathlib;
END exp; END exp;
PROCEDURE longexp(x: LONGREAL): LONGREAL; PROCEDURE longexp(x: LONGREAL): LONGREAL;
(* (* 2**x = (Q(x*x)+x*P(x*x))/(Q(x*x)-x*P(x*x)) for x in [0,0.5] *)
* n = floor(x / ln2), d = x / ln2 - n (* Hart & Cheney #1069 *)
* exp(x) = exp((x / ln2) * ln2) = exp((n + d) * ln2) =
* exp(n * ln2) * exp(d * ln2) = 2 ** n * exp(d * ln2)
*)
CONST CONST
a1 = -0.9999999995D; p0 = 0.2080384346694663001443843411D+07;
a2 = 0.4999999206D; p1 = 0.3028697169744036299076048876D+05;
a3 = -0.1666653019D; p2 = 0.6061485330061080841615584556D+02;
a4 = 0.0416573475D; q0 = 0.6002720360238832528230907598D+07;
a5 = -0.0083013598D; q1 = 0.3277251518082914423057964422D+06;
a6 = 0.0013298820D; q2 = 0.1749287689093076403844945335D+04;
a7 = -0.0001413161D; q3 = 0.1000000000000000000000000000D+01;
VAR VAR
neg: BOOLEAN; neg: BOOLEAN;
polval: LONGREAL; xPxx, Qxx: LONGREAL;
n: LONGREAL; n: LONGREAL;
n1 : INTEGER; n1 : INTEGER;
xsq : LONGREAL;
large: BOOLEAN;
BEGIN BEGIN
neg := x < 0.0D; neg := x < 0.0D;
IF neg THEN IF neg THEN
x := -x; x := -x;
END; END;
x := FIF(x, 1.0D/LONG(ln2), n) * LONG(ln2); x := FIF(x/longln2, 1.0D, n);
polval := 1.0D /(1.0D + x*(a1+x*(a2+x*(a3+x*(a4+x*(a5+x*(a6+x*a7))))))); large := x > 0.5D;
IF large THEN x := x - 0.5D; END;
xsq := x*x;
xPxx := x*((p2*xsq+p1)*xsq+p0);
Qxx := ((q3*xsq+q2)*xsq+q1)*xsq+q0;
x := (Qxx + xPxx)/(Qxx - xPxx);
IF large THEN
x := x * Sqrt2;
END;
n1 := TRUNCD(n + 0.5D); n1 := TRUNCD(n + 0.5D);
WHILE n1 >= 16 DO WHILE n1 >= 16 DO
polval := polval * 65536.0D; x := x * 65536.0D;
n1 := n1 - 16; n1 := n1 - 16;
END; END;
WHILE n1 > 0 DO WHILE n1 > 0 DO
polval := polval * 2.0D; x := x * 2.0D;
DEC(n1); DEC(n1);
END; END;
IF neg THEN RETURN 1.0D/polval; END; IF neg THEN RETURN 1.0D/x; END;
RETURN polval; RETURN x;
END longexp; END longexp;
PROCEDURE ln(x: REAL): REAL; (* natural log *) PROCEDURE ln(x: REAL): REAL; (* natural log *)
@ -149,18 +152,23 @@ IMPLEMENTATION MODULE Mathlib;
END ln; END ln;
PROCEDURE longln(x: LONGREAL): LONGREAL; (* natural log *) PROCEDURE longln(x: LONGREAL): LONGREAL; (* natural log *)
(* log(x) = z*P(z*z)/Q(z*z), z = (x-1)/(x+1), x in [1/sqrt(2), sqrt(2)]
Hart & Cheney #2707
*)
CONST CONST
a1 = 0.9999964239D; p0 = 0.7504094990777122217455611007D+02;
a2 = -0.4998741238D; p1 = -0.1345669115050430235318253537D+03;
a3 = 0.3317990258D; p2 = 0.7413719213248602512779336470D+02;
a4 = -0.2407338084D; p3 = -0.1277249755012330819984385000D+02;
a5 = 0.1676540711D; p4 = 0.3327108381087686938144000000D+00;
a6 = -0.0953293897D; q0 = 0.3752047495388561108727775374D+02;
a7 = 0.0360884937D; q1 = -0.7979028073715004879439951583D+02;
a8 = -0.0064535442D; q2 = 0.5616126132118257292058560360D+02;
q3 = -0.1450868091858082685362325000D+02;
q4 = 0.1000000000000000000000000000D+01;
VAR VAR
exp: INTEGER; exp: INTEGER;
polval: LONGREAL; z, zsq: LONGREAL;
BEGIN BEGIN
IF x <= 0.0D THEN IF x <= 0.0D THEN
@ -168,13 +176,15 @@ IMPLEMENTATION MODULE Mathlib;
HALT HALT
END; END;
x := FEF(x, exp); x := FEF(x, exp);
WHILE x < 1.0D DO WHILE x < OneOverSqrt2 DO
x := x + x; x := x + x;
DEC(exp); DEC(exp);
END; END;
x := x - 1.0D; z := (x - 1.0D) / (x + 1.0D);
polval := x*(a1+x*(a2+x*(a3+x*(a4+x*(a5+x*(a6+x*(a7+a8*x))))))); zsq := z*z;
RETURN polval + FLOATD(exp) * LONG(ln2); RETURN z * ((((p4*zsq+p3)*zsq+p2)*zsq+p1)*zsq+p0) /
((((q4*zsq+q3)*zsq+q2)*zsq+q1)*zsq+q0) +
FLOATD(exp) * longln2;
END longln; END longln;
PROCEDURE log(x: REAL): REAL; (* log with base 10 *) PROCEDURE log(x: REAL): REAL; (* log with base 10 *)
@ -184,7 +194,7 @@ IMPLEMENTATION MODULE Mathlib;
PROCEDURE longlog(x: LONGREAL): LONGREAL; (* log with base 10 *) PROCEDURE longlog(x: LONGREAL): LONGREAL; (* log with base 10 *)
BEGIN BEGIN
RETURN longln(x)/LONG(ln10); RETURN longln(x)/longln10;
END longlog; END longlog;
(* trigonometric functions; arguments in radians *) (* trigonometric functions; arguments in radians *)
@ -194,34 +204,80 @@ IMPLEMENTATION MODULE Mathlib;
RETURN SHORT(longsin(LONG(x))); RETURN SHORT(longsin(LONG(x)));
END sin; END sin;
PROCEDURE longsin(x: LONGREAL): LONGREAL; PROCEDURE sinus(x: LONGREAL; quadrant: INTEGER) : LONGREAL;
(* sin(0.5*pi*x) = x * P(x*x)/Q(x*x) for x in [0,1]
Hart & Cheney # 3374
*)
CONST CONST
a2 = -0.1666666664D; p0 = 0.4857791909822798473837058825D+10;
a4 = 0.0083333315D; p1 = -0.1808816670894030772075877725D+10;
a6 = -0.0001984090D; p2 = 0.1724314784722489597789244188D+09;
a8 = 0.0000027526D; p3 = -0.6351331748520454245913645971D+07;
a10 = -0.0000000239D; p4 = 0.1002087631419532326179108883D+06;
p5 = -0.5830988897678192576148973679D+03;
q0 = 0.3092566379840468199410228418D+10;
q1 = 0.1202384907680254190870913060D+09;
q2 = 0.2321427631602460953669856368D+07;
q3 = 0.2848331644063908832127222835D+05;
q4 = 0.2287602116741682420054505174D+03;
q5 = 0.1000000000000000000000000000D+01;
A1 = 6.2822265625D;
A2 = 0.00095874467958647692528676655900576D;
VAR VAR
xsqr: LONGREAL; xsq, x1, x2, n : LONGREAL;
neg: BOOLEAN; t : INTEGER;
BEGIN BEGIN
neg := FALSE;
IF x < 0.0D THEN IF x < 0.0D THEN
neg := TRUE; INC(quadrant, 2);
x := -x; x := -x;
END; END;
x := FIF(x, 1.0D / LONG(twicepi), (* dummy *) xsqr) * LONG(twicepi); IF longhalfpi - x = longhalfpi THEN
IF x >= LONG(pi) THEN CASE quadrant OF
neg := NOT neg; | 0,2:
x := x - LONG(pi); RETURN 0.0D;
| 1:
RETURN 1.0D;
| 3:
RETURN -1.0D;
END;
END; END;
IF x > LONG(halfpi) THEN IF x >= longtwicepi THEN
x := LONG(pi) - x; IF x <= FLOATD(MAX(LONGINT)) THEN
(* Use extended precision to calculate reduced argument.
Split 2pi in 2 parts a1 and a2, of which the first only
uses some bits of the mantissa, so that n * a1 is
exactly representable, where n is the integer part of
x/pi.
Here we used 12 bits of the mantissa for a1.
Also split x in integer part x1 and fraction part x2.
We then compute x-n*2pi as ((x1 - n*a1) + x2) - n*a2.
*)
n := FLOATD(TRUNCD(x/longtwicepi));
x1 := FLOATD(TRUNCD(x));
x2 := x - x1;
x := ((x1 - n * A1) + x2) - n * A2;
ELSE
x := FIF(x/longtwicepi, 1.0D, x1) * longtwicepi;
END
END; END;
xsqr := x * x; x := x / longhalfpi;
x := x * (1.0D + xsqr*(a2+xsqr*(a4+xsqr*(a6+xsqr*(a8+xsqr*a10))))); t := TRUNC(x);
IF neg THEN RETURN -x; END; x := x - FLOATD(t);
RETURN x; quadrant := (quadrant + t MOD 4) MOD 4;
IF ODD(quadrant) THEN
x := 1.0D - x;
END;
IF quadrant > 1 THEN
x := -x;
END;
xsq := x * x;
RETURN x * (((((p5*xsq+p4)*xsq+p3)*xsq+p2)*xsq+p1)*xsq+p0) /
(((((q5*xsq+q4)*xsq+q3)*xsq+q2)*xsq+q1)*xsq+q0);
END sinus;
PROCEDURE longsin(x: LONGREAL): LONGREAL;
BEGIN
RETURN sinus(x, 0);
END longsin; END longsin;
PROCEDURE cos(x: REAL): REAL; PROCEDURE cos(x: REAL): REAL;
@ -230,30 +286,9 @@ IMPLEMENTATION MODULE Mathlib;
END cos; END cos;
PROCEDURE longcos(x: LONGREAL): LONGREAL; PROCEDURE longcos(x: LONGREAL): LONGREAL;
CONST
a2 = -0.4999999963D;
a4 = 0.0416666418D;
a6 = -0.0013888397D;
a8 = 0.0000247609D;
a10 = -0.0000002605D;
VAR
xsqr: LONGREAL;
neg: BOOLEAN;
BEGIN BEGIN
neg := FALSE;
IF x < 0.0D THEN x := -x; END; IF x < 0.0D THEN x := -x; END;
x := FIF(x, 1.0D / LONG(twicepi), (* dummy *) xsqr) * LONG(twicepi); RETURN sinus(x, 1);
IF x >= LONG(pi) THEN
x := LONG(twicepi) - x;
END;
IF x > LONG(halfpi) THEN
neg := NOT neg;
x := LONG(pi) - x;
END;
xsqr := x * x;
x := 1.0D + xsqr*(a2+xsqr*(a4+xsqr*(a6+xsqr*(a8+xsqr*a10))));
IF neg THEN RETURN -x; END;
RETURN x;
END longcos; END longcos;
PROCEDURE tan(x: REAL): REAL; PROCEDURE tan(x: REAL): REAL;
@ -277,24 +312,31 @@ IMPLEMENTATION MODULE Mathlib;
RETURN SHORT(longarcsin(LONG(x))); RETURN SHORT(longarcsin(LONG(x)));
END arcsin; END arcsin;
PROCEDURE longarcsin(x: LONGREAL): LONGREAL; PROCEDURE arcsincos(x: LONGREAL; cosfl: BOOLEAN): LONGREAL;
CONST VAR
a0 = 1.5707963050D; negative : BOOLEAN;
a1 = -0.2145988016D;
a2 = 0.0889789874D;
a3 = -0.0501743046D;
a4 = 0.0308918810D;
a5 = -0.0170881256D;
a6 = 0.0066700901D;
a7 = -0.0012624911D;
BEGIN BEGIN
IF x < 0.0D THEN x := -x; END; negative := x <= 0.0D;
IF negative THEN x := -x; END;
IF x > 1.0D THEN IF x > 1.0D THEN
Message("arcsin: argument > 1"); Message("arcsin or arccos: argument > 1");
HALT HALT
END; END;
RETURN LONG(halfpi) - IF x = 1.0D THEN
longsqrt(1.0D - x)*(a0+x*(a1+x*(a2+x*(a3+x*(a4+x*(a5+x*(a6+x*a7))))))); x := longhalfpi;
ELSE
x := longarctan(x/longsqrt(1.0D - x*x));
END;
IF negative THEN x := -x; END;
IF cosfl THEN
RETURN longhalfpi - x;
END;
RETURN x;
END arcsincos;
PROCEDURE longarcsin(x: LONGREAL): LONGREAL;
BEGIN
RETURN arcsincos(x, FALSE);
END longarcsin; END longarcsin;
PROCEDURE arccos(x: REAL): REAL; PROCEDURE arccos(x: REAL): REAL;
@ -304,7 +346,7 @@ IMPLEMENTATION MODULE Mathlib;
PROCEDURE longarccos(x: LONGREAL): LONGREAL; PROCEDURE longarccos(x: LONGREAL): LONGREAL;
BEGIN BEGIN
RETURN LONG(halfpi) - longarcsin(x); RETURN arcsincos(x, TRUE);
END longarccos; END longarccos;
PROCEDURE arctan(x: REAL): REAL; PROCEDURE arctan(x: REAL): REAL;
@ -312,35 +354,109 @@ IMPLEMENTATION MODULE Mathlib;
RETURN SHORT(longarctan(LONG(x))); RETURN SHORT(longarctan(LONG(x)));
END arctan; END arctan;
TYPE
precomputed = RECORD
X: LONGREAL; (* partition point *)
arctan: LONGREAL; (* arctan of evaluation node *)
OneOverXn: LONGREAL; (* 1/xn *)
OneOverXnSquarePlusone: LONGREAL; (* ... *)
END;
VAR arctaninit: BOOLEAN;
precomp : ARRAY[0..4] OF precomputed;
PROCEDURE longarctan(x: LONGREAL): LONGREAL; PROCEDURE longarctan(x: LONGREAL): LONGREAL;
(* The interval [0, infinity) is treated as follows:
Define partition points Xi
X0 = 0
X1 = tan(pi/16)
X2 = tan(3pi/16)
X3 = tan(5pi/16)
X4 = tan(7pi/16)
X5 = infinity
and evaluation nodes xi
x2 = tan(2pi/16)
x3 = tan(4pi/16)
x4 = tan(6pi/16)
x5 = infinity
An argument x in [Xn-1, Xn] is now reduced to an argument
t in [-X1, X1] by the following formulas:
t = 1/xn - (1/(xn*xn) + 1)/((1/xn) + x)
arctan(x) = arctan(xi) + arctan(t)
For the interval [0, p/16] an approximation is used:
arctan(x) = x * P(x*x)/Q(x*x)
*)
(* Hart & Cheney # 5037 *)
CONST CONST
a2 = -0.3333314528D; p0 = 0.7698297257888171026986294745D+03;
a4 = 0.1999355085D; p1 = 0.1557282793158363491416585283D+04;
a6 = -0.1420889944D; p2 = 0.1033384651675161628243434662D+04;
a8 = 0.1065626393D; p3 = 0.2485841954911840502660889866D+03;
a10 = -0.0752896400D; p4 = 0.1566564964979791769948970100D+02;
a12 = 0.0429096318D; q0 = 0.7698297257888171026986294911D+03;
a14 = -0.0161657367D; q1 = 0.1813892701754635858982709369D+04;
a16 = 0.0028662257D; q2 = 0.1484049607102276827437401170D+04;
q3 = 0.4904645326203706217748848797D+03;
q4 = 0.5593479839280348664778328000D+02;
q5 = 0.1000000000000000000000000000D+01;
VAR VAR
xsqr: LONGREAL; xsqr: LONGREAL;
rev: BOOLEAN;
neg: BOOLEAN; neg: BOOLEAN;
i: INTEGER;
BEGIN BEGIN
rev := FALSE; IF NOT arctaninit THEN
arctaninit := TRUE;
WITH precomp[0] DO
X := 0.19891236737965800691159762264467622;
arctan := 0.0D;
OneOverXn := 0.0D;
OneOverXnSquarePlusone := 0.0D;
END;
WITH precomp[1] DO
X := 0.66817863791929891999775768652308076;
arctan := longpi/8.0D;
OneOverXn := 2.41421356237309504880168872420969808;
OneOverXnSquarePlusone := 6.82842712474619009760337744841939616;
END;
WITH precomp[2] DO
X := 1.49660576266548901760113513494247691;
arctan := longquartpi;
OneOverXn := 1.0;
OneOverXnSquarePlusone := 2.0;
END;
WITH precomp[3] DO
X := 5.02733949212584810451497507106407238;
arctan := 3.0D*longpi/8.0D;
OneOverXn := 0.41421356237309504880168872420969808;
OneOverXnSquarePlusone := 1.17157287525380998659662255158060384;
END;
WITH precomp[4] DO
X := 0.0;
arctan := longhalfpi;
OneOverXn := 0.0;
OneOverXnSquarePlusone := 1.0;
END;
END;
neg := FALSE; neg := FALSE;
IF x < 0.0D THEN IF x < 0.0D THEN
neg := TRUE; neg := TRUE;
x := -x; x := -x;
END; END;
IF x > 1.0D THEN i := 0;
rev := TRUE; WHILE (i <= 3) AND (x <= precomp[i].X) DO
x := 1.0D / x; INC(i);
END; END;
xsqr := x * x; IF (i # 0) THEN
x := x * (1.0D + xsqr*(a2+xsqr*(a4+xsqr*(a6+xsqr*(a8+xsqr*(a10+xsqr*(a12+xsqr*(a14+xsqr*a16)))))))); WITH precomp[i] DO
IF rev THEN x := arctan + longarctan(OneOverXn-OneOverXnSquarePlusone/(OneOverXn+x));
x := LONG(quartpi) - x; END
ELSE
xsqr := x * x;
x := x * ((((p4*xsqr+p3)*xsqr+p2)*xsqr+p1)*xsqr+p0) /
(((((q5*xsqr+q4)*xsqr+q3)*xsqr+q2)*xsqr+q1)*xsqr+q0);
END; END;
IF neg THEN RETURN -x; END; IF neg THEN RETURN -x; END;
RETURN x; RETURN x;
@ -452,4 +568,6 @@ IMPLEMENTATION MODULE Mathlib;
RETURN x * OneDegreeInRadians; RETURN x * OneDegreeInRadians;
END longDegreeToRadian; END longDegreeToRadian;
BEGIN
arctaninit := FALSE;
END Mathlib. END Mathlib.