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

use better algorithms for some mathematical functions

Browse source
This commit is contained in:
ceriel 1989-06-20 13:10:32 +00:00
parent 6d78cd6710
commit 11349c78cd
1 changed files with 229 additions and 226 deletions

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

@ -17,7 +17,6 @@ IMPLEMENTATION MODULE Mathlib;
CONST CONST
OneRadianInDegrees = 57.295779513082320876798155D; OneRadianInDegrees = 57.295779513082320876798155D;
OneDegreeInRadians = 0.017453292519943295769237D; OneDegreeInRadians = 0.017453292519943295769237D;
Sqrt2 = 1.41421356237309504880168872420969808D;
OneOverSqrt2 = 0.70710678118654752440084436210484904D; OneOverSqrt2 = 0.70710678118654752440084436210484904D;
(* basic functions *) (* basic functions *)
@ -94,56 +93,68 @@ IMPLEMENTATION MODULE Mathlib;
RETURN temp; RETURN temp;
END longsqrt; END longsqrt;
PROCEDURE ldexp(x:LONGREAL; n: INTEGER): LONGREAL;
BEGIN
WHILE n >= 16 DO
x := x * 65536.0D;
n := n - 16;
END;
WHILE n > 0 DO
x := x * 2.0D;
DEC(n);
END;
WHILE n <= -16 DO
x := x / 65536.0D;
n := n + 16;
END;
WHILE n < 0 DO
x := x / 2.0D;
INC(n);
END;
RETURN x;
END ldexp;
PROCEDURE exp(x: REAL): REAL; PROCEDURE exp(x: REAL): REAL;
BEGIN BEGIN
RETURN SHORT(longexp(LONG(x))); RETURN SHORT(longexp(LONG(x)));
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] *) (* Algorithm and coefficients from:
(* Hart & Cheney #1069 *) "Software manual for the elementary functions"
by W.J. Cody and W. Waite, Prentice-Hall, 1980
*)
CONST CONST
p0 = 0.2080384346694663001443843411D+07; p0 = 0.25000000000000000000D+00;
p1 = 0.3028697169744036299076048876D+05; p1 = 0.75753180159422776666D-02;
p2 = 0.6061485330061080841615584556D+02; p2 = 0.31555192765684646356D-04;
q0 = 0.6002720360238832528230907598D+07; q0 = 0.50000000000000000000D+00;
q1 = 0.3277251518082914423057964422D+06; q1 = 0.56817302698551221787D-01;
q2 = 0.1749287689093076403844945335D+04; q2 = 0.63121894374398503557D-03;
q3 = 0.1000000000000000000000000000D+01; q3 = 0.75104028399870046114D-06;
VAR VAR
neg: BOOLEAN; neg: BOOLEAN;
xPxx, Qxx: LONGREAL; n: INTEGER;
n: LONGREAL; xn, g, x1, x2: LONGREAL;
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/longln2, 1.0D, n); n := TRUNC(x/longln2 + 0.5D);
large := x > 0.5D; xn := FLOATD(n);
IF large THEN x := x - 0.5D; END; x1 := FLOATD(TRUNCD(x));
xsq := x*x; x2 := x - x1;
xPxx := x*((p2*xsq+p1)*xsq+p0); g := ((x1 - xn * 0.693359375D)+x2) - xn * (-2.1219444005469058277D-4);
Qxx := ((q3*xsq+q2)*xsq+q1)*xsq+q0; IF neg THEN
x := (Qxx + xPxx)/(Qxx - xPxx); g := -g;
IF large THEN n := -n;
x := x * Sqrt2;
END; END;
n1 := TRUNCD(n + 0.5D); xn := g*g;
WHILE n1 >= 16 DO x := g*((p2*xn+p1)*xn+p0);
x := x * 65536.0D; INC(n);
n1 := n1 - 16; RETURN ldexp(0.5D + x/((((q3*xn+q2)*xn+q1)*xn+q0) - x), n);
END;
WHILE n1 > 0 DO
x := x * 2.0D;
DEC(n1);
END;
IF neg THEN RETURN 1.0D/x; END;
RETURN x;
END longexp; END longexp;
PROCEDURE ln(x: REAL): REAL; (* natural log *) PROCEDURE ln(x: REAL): REAL; (* natural log *)
@ -152,23 +163,21 @@ 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)] (* Algorithm and coefficients from:
Hart & Cheney #2707 "Software manual for the elementary functions"
*) by W.J. Cody and W. Waite, Prentice-Hall, 1980
*)
CONST CONST
p0 = 0.7504094990777122217455611007D+02; p0 = -0.64124943423745581147D+02;
p1 = -0.1345669115050430235318253537D+03; p1 = 0.16383943563021534222D+02;
p2 = 0.7413719213248602512779336470D+02; p2 = -0.78956112887491257267D+00;
p3 = -0.1277249755012330819984385000D+02; q0 = -0.76949932108494879777D+03;
p4 = 0.3327108381087686938144000000D+00; q1 = 0.31203222091924532844D+03;
q0 = 0.3752047495388561108727775374D+02; q2 = -0.35667977739034646171D+02;
q1 = -0.7979028073715004879439951583D+02; q3 = 1.0D;
q2 = 0.5616126132118257292058560360D+02;
q3 = -0.1450868091858082685362325000D+02;
q4 = 0.1000000000000000000000000000D+01;
VAR VAR
exp: INTEGER; exp: INTEGER;
z, zsq: LONGREAL; z, znum, zden, w: LONGREAL;
BEGIN BEGIN
IF x <= 0.0D THEN IF x <= 0.0D THEN
@ -176,15 +185,20 @@ IMPLEMENTATION MODULE Mathlib;
HALT HALT
END; END;
x := FEF(x, exp); x := FEF(x, exp);
WHILE x < OneOverSqrt2 DO IF x > OneOverSqrt2 THEN
x := x + x; znum := (x - 0.5D) - 0.5D;
zden := x * 0.5D + 0.5D;
ELSE
znum := x - 0.5D;
zden := znum * 0.5D + 0.5D;
DEC(exp); DEC(exp);
END; END;
z := (x - 1.0D) / (x + 1.0D); z := znum / zden;
zsq := z*z; w := z * z;
RETURN z * ((((p4*zsq+p3)*zsq+p2)*zsq+p1)*zsq+p0) / x := z + z * w * (((p2*w+p1)*w+p0)/(((q3*w+q2)*w+q1)*w+q0));
((((q4*zsq+q3)*zsq+q2)*zsq+q1)*zsq+q0) + z := FLOATD(exp);
FLOATD(exp) * longln2; x := x + z * (-2.121944400546905827679D-4);
RETURN x + z * 0.693359375D;
END longln; END longln;
PROCEDURE log(x: REAL): REAL; (* log with base 10 *) PROCEDURE log(x: REAL): REAL; (* log with base 10 *)
@ -204,80 +218,60 @@ IMPLEMENTATION MODULE Mathlib;
RETURN SHORT(longsin(LONG(x))); RETURN SHORT(longsin(LONG(x)));
END sin; END sin;
PROCEDURE sinus(x: LONGREAL; quadrant: INTEGER) : LONGREAL; PROCEDURE sinus(x: LONGREAL; cosflag: BOOLEAN) : LONGREAL;
(* sin(0.5*pi*x) = x * P(x*x)/Q(x*x) for x in [0,1] (* Algorithm and coefficients from:
Hart & Cheney # 3374 "Software manual for the elementary functions"
by W.J. Cody and W. Waite, Prentice-Hall, 1980
*) *)
CONST CONST
p0 = 0.4857791909822798473837058825D+10; r0 = -0.16666666666666665052D+00;
p1 = -0.1808816670894030772075877725D+10; r1 = 0.83333333333331650314D-02;
p2 = 0.1724314784722489597789244188D+09; r2 = -0.19841269841201840457D-03;
p3 = -0.6351331748520454245913645971D+07; r3 = 0.27557319210152756119D-05;
p4 = 0.1002087631419532326179108883D+06; r4 = -0.25052106798274584544D-07;
p5 = -0.5830988897678192576148973679D+03; r5 = 0.16058936490371589114D-09;
q0 = 0.3092566379840468199410228418D+10; r6 = -0.76429178068910467734D-12;
q1 = 0.1202384907680254190870913060D+09; r7 = 0.27204790957888846175D-14;
q2 = 0.2321427631602460953669856368D+07; A1 = 3.1416015625D;
q3 = 0.2848331644063908832127222835D+05; A2 = -8.908910206761537356617D-6;
q4 = 0.2287602116741682420054505174D+03;
q5 = 0.1000000000000000000000000000D+01;
A1 = 6.2822265625D;
A2 = 0.00095874467958647692528676655900576D;
VAR VAR
xsq, x1, x2, n : LONGREAL; x1, x2, y : LONGREAL;
t : INTEGER; neg : BOOLEAN;
BEGIN BEGIN
IF x < 0.0D THEN IF x < 0.0D THEN
INC(quadrant, 2); neg := TRUE;
x := -x; x := -x
ELSE neg := FALSE
END; END;
IF longhalfpi - x = longhalfpi THEN IF cosflag THEN
CASE quadrant OF neg := FALSE;
| 0,2: y := longhalfpi + x
RETURN 0.0D; ELSE
| 1: y := x
RETURN 1.0D;
| 3:
RETURN -1.0D;
END;
END; END;
IF x >= longtwicepi THEN y := y / longpi + 0.5D;
IF x <= FLOATD(MAX(LONGINT)) THEN
(* Use extended precision to calculate reduced argument. IF FIF(y, 1.0D, y) < 0.0D THEN ; END;
Split 2pi in 2 parts a1 and a2, of which the first only IF FIF(y, 0.5D, x1) # 0.0D THEN neg := NOT neg END;
uses some bits of the mantissa, so that n * a1 is IF cosflag THEN y := y - 0.5D END;
exactly representable, where n is the integer part of x2 := FIF(x, 1.0, x1);
x/pi. x := x1 - y * A1;
Here we used 12 bits of the mantissa for a1. x := x + x2;
Also split x in integer part x1 and fraction part x2. x := x - y * A2;
We then compute x-n*2pi as ((x1 - n*a1) + x2) - n*a2.
*) IF x < 0.0D THEN
n := FLOATD(TRUNCD(x/longtwicepi)); neg := NOT neg;
x1 := FLOATD(TRUNCD(x)); x := -x
x2 := x - x1;
x := ((x1 - n * A1) + x2) - n * A2;
ELSE
x := FIF(x/longtwicepi, 1.0D, x1) * longtwicepi;
END
END; END;
x := x / longhalfpi; y := x * x;
t := TRUNC(x); x := x + x * y * (((((((r7*y+r6)*y+r5)*y+r4)*y+r3)*y+r2)*y+r1)*y+r0);
x := x - FLOATD(t); IF neg THEN RETURN -x END;
quadrant := (quadrant + t MOD 4) MOD 4; RETURN x;
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; END sinus;
PROCEDURE longsin(x: LONGREAL): LONGREAL; PROCEDURE longsin(x: LONGREAL): LONGREAL;
BEGIN BEGIN
RETURN sinus(x, 0); RETURN sinus(x, FALSE);
END longsin; END longsin;
PROCEDURE cos(x: REAL): REAL; PROCEDURE cos(x: REAL): REAL;
@ -288,7 +282,7 @@ IMPLEMENTATION MODULE Mathlib;
PROCEDURE longcos(x: LONGREAL): LONGREAL; PROCEDURE longcos(x: LONGREAL): LONGREAL;
BEGIN BEGIN
IF x < 0.0D THEN x := -x; END; IF x < 0.0D THEN x := -x; END;
RETURN sinus(x, 1); RETURN sinus(x, TRUE);
END longcos; END longcos;
PROCEDURE tan(x: REAL): REAL; PROCEDURE tan(x: REAL): REAL;
@ -297,14 +291,49 @@ IMPLEMENTATION MODULE Mathlib;
END tan; END tan;
PROCEDURE longtan(x: LONGREAL): LONGREAL; PROCEDURE longtan(x: LONGREAL): LONGREAL;
VAR cosinus: LONGREAL; (* Algorithm and coefficients from:
"Software manual for the elementary functions"
by W.J. Cody and W. Waite, Prentice-Hall, 1980
*)
CONST
p1 = -0.13338350006421960681D+00;
p2 = 0.34248878235890589960D-02;
p3 = -0.17861707342254426711D-04;
q0 = 1.0D;
q1 = -0.46671683339755294240D+00;
q2 = 0.25663832289440112864D-01;
q3 = -0.31181531907010027307D-03;
q4 = 0.49819433993786512270D-06;
A1 = 1.57080078125D;
A2 = -4.454455103380768678308D-06;
VAR y, x1, x2: LONGREAL;
negative: BOOLEAN;
invert: BOOLEAN;
BEGIN BEGIN
cosinus := longcos(x); negative := x < 0.0D;
IF cosinus = 0.0D THEN y := x / longhalfpi + 0.5D;
Message("tan: result does not exist");
HALT (* Use extended precision to calculate reduced argument.
END; Here we used 12 bits of the mantissa for a1.
RETURN longsin(x)/cosinus; Also split x in integer part x1 and fraction part x2.
*)
IF FIF(y, 1.0D, y) < 0.0D THEN ; END;
invert := FIF(y, 0.5D, x1) # 0.0D;
x2 := FIF(x, 1.0D, x1);
x := x1 - y * A1;
x := x + x2;
x := x - y * A2;
y := x * x;
x := x + x * y * ((p3*y+p2)*y+p1);
y := (((q4*y+q3)*y+q2)*y+q1)*y+q0;
IF negative THEN x := -x END;
IF invert THEN RETURN -y/x END;
RETURN x/y;
END longtan; END longtan;
PROCEDURE arcsin(x: REAL): REAL; PROCEDURE arcsin(x: REAL): REAL;
@ -313,24 +342,48 @@ IMPLEMENTATION MODULE Mathlib;
END arcsin; END arcsin;
PROCEDURE arcsincos(x: LONGREAL; cosfl: BOOLEAN): LONGREAL; PROCEDURE arcsincos(x: LONGREAL; cosfl: BOOLEAN): LONGREAL;
VAR CONST
p0 = -0.27368494524164255994D+02;
p1 = 0.57208227877891731407D+02;
p2 = -0.39688862997540877339D+02;
p3 = 0.10152522233806463645D+02;
p4 = -0.69674573447350646411D+00;
q0 = -0.16421096714498560795D+03;
q1 = 0.41714430248260412556D+03;
q2 = -0.38186303361750149284D+03;
q3 = 0.15095270841030604719D+03;
q4 = -0.23823859153670238830D+02;
q5 = 1.0D;
VAR
negative : BOOLEAN; negative : BOOLEAN;
big: BOOLEAN;
g: LONGREAL;
BEGIN BEGIN
negative := x <= 0.0D; negative := x < 0.0D;
IF negative THEN x := -x; END; IF negative THEN x := -x; END;
IF x > 1.0D THEN IF x > 0.5D THEN
Message("arcsin or arccos: argument > 1"); big := TRUE;
HALT IF x > 1.0D THEN
END; Message("arcsin or arccos: argument > 1");
IF x = 1.0D THEN HALT
x := longhalfpi; END;
g := 0.5D - 0.5D * x;
x := -longsqrt(g);
x := x + x;
ELSE ELSE
x := longarctan(x/longsqrt(1.0D - x*x)); big := FALSE;
g := x * x;
END; END;
IF negative THEN x := -x; END; x := x + x * g *
IF cosfl THEN ((((p4*g+p3)*g+p2)*g+p1)*g+p0)/(((((q5*g+q4)*g+q3)*g+q2)*g+q1)*g+q0);
RETURN longhalfpi - x; IF cosfl AND NOT negative THEN x := -x END;
IF cosfl = NOT big THEN
x := (x + longquartpi) + longquartpi;
ELSIF cosfl AND negative AND big THEN
x := (x + longhalfpi) + longhalfpi;
END; END;
IF negative AND NOT cosfl THEN x := -x END;
RETURN x; RETURN x;
END arcsincos; END arcsincos;
@ -354,115 +407,65 @@ IMPLEMENTATION MODULE Mathlib;
RETURN SHORT(longarctan(LONG(x))); RETURN SHORT(longarctan(LONG(x)));
END arctan; END arctan;
TYPE VAR A: ARRAY[0..3] OF LONGREAL;
precomputed = RECORD arctaninit: BOOLEAN;
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: (* Algorithm and coefficients from:
Define partition points Xi "Software manual for the elementary functions"
X0 = 0 by W.J. Cody and W. Waite, Prentice-Hall, 1980
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, tan(p/16)] an approximation is used:
arctan(x) = x * P(x*x)/Q(x*x)
*) *)
(* Hart & Cheney # 5037 *)
CONST CONST
p0 = 0.7698297257888171026986294745D+03; p0 = -0.13688768894191926929D+02;
p1 = 0.1557282793158363491416585283D+04; p1 = -0.20505855195861651981D+02;
p2 = 0.1033384651675161628243434662D+04; p2 = -0.84946240351320683534D+01;
p3 = 0.2485841954911840502660889866D+03; p3 = -0.83758299368150059274D+00;
p4 = 0.1566564964979791769948970100D+02; q0 = 0.41066306682575781263D+02;
q0 = 0.7698297257888171026986294911D+03; q1 = 0.86157349597130242515D+02;
q1 = 0.1813892701754635858982709369D+04; q2 = 0.59578436142597344465D+02;
q2 = 0.1484049607102276827437401170D+04; q3 = 0.15024001160028576121D+02;
q3 = 0.4904645326203706217748848797D+03; q4 = 1.0D;
q4 = 0.5593479839280348664778328000D+02;
q5 = 0.1000000000000000000000000000D+01;
VAR VAR
xsqr: LONGREAL; g: LONGREAL;
neg: BOOLEAN; neg: BOOLEAN;
i: INTEGER; n: INTEGER;
BEGIN BEGIN
IF NOT arctaninit THEN IF NOT arctaninit THEN
arctaninit := TRUE; arctaninit := TRUE;
WITH precomp[0] DO A[0] := 0.0D;
X := 0.19891236737965800691159762264467622D; A[1] := 0.52359877559829887307710723554658381D; (* p1/6 *)
arctan := 0.0D; A[2] := longhalfpi;
OneOverXn := 0.0D; A[3] := 1.04719755119659774615421446109316763D; (* pi/3 *)
OneOverXnSquarePlusone := 0.0D;
END;
WITH precomp[1] DO
X := 0.66817863791929891999775768652308076D;
arctan := 0.39269908169872415480783042290993786D;
OneOverXn := 2.41421356237309504880168872420969808D;
OneOverXnSquarePlusone := 6.82842712474619009760337744841939616D;
END;
WITH precomp[2] DO
X := 1.49660576266548901760113513494247691D;
arctan := longquartpi;
OneOverXn := 1.0;
OneOverXnSquarePlusone := 2.0;
END;
WITH precomp[3] DO
X := 5.02733949212584810451497507106407238D;
arctan := 1.17809724509617246442349126872981358D;
OneOverXn := 0.41421356237309504880168872420969808D;
OneOverXnSquarePlusone := 1.17157287525380998659662255158060384D;
END;
WITH precomp[4] DO
X := 0.0D;
arctan := longhalfpi;
OneOverXn := 0.0D;
OneOverXnSquarePlusone := 1.0D;
END;
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;
i := 0; IF x > 1.0D THEN
WHILE (i <= 3) AND (x >= precomp[i].X) DO x := 1.0D/x;
INC(i); n := 2
END;
IF (i # 0) THEN
WITH precomp[i] DO
x := arctan + longarctan(OneOverXn-OneOverXnSquarePlusone/(OneOverXn+x));
END
ELSE ELSE
xsqr := x * x; n := 0
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 x > 0.26794919243112270647D (* 2-sqrt(3) *) THEN
INC(n);
x := (((0.73205080756887729353D*x-0.5D)-0.5D)+x)/
(1.73205080756887729353D + x);
END;
g := x*x;
x := x + x * g * (((p3*g+p2)*g+p1)*g+p0) / ((((q4*g+q3)*g+q2)*g+q1)*g+q0);
IF n > 1 THEN x := -x END;
x := x + A[n];
IF neg THEN RETURN -x; END; IF neg THEN RETURN -x; END;
RETURN x; RETURN x;
END longarctan; END longarctan;
(* hyperbolic functions *) (* hyperbolic functions *)
(* The C math library has better implementations for some of these, but
they depend on some properties of the floating point implementation,
and, for now, we don't want that in the Modula-2 system.
*)
PROCEDURE sinh(x: REAL): REAL; PROCEDURE sinh(x: REAL): REAL;
BEGIN BEGIN