From dbbff76f4c12776d335a6b84a0dfc77543804e4d Mon Sep 17 00:00:00 2001 From: ceriel Date: Mon, 25 Jul 1988 16:41:51 +0000 Subject: [PATCH] Used new math lib of C to create new version of Mathlib --- lang/m2/libm2/Mathlib.def | 30 ++-- lang/m2/libm2/Mathlib.mod | 352 +++++++++++++++++++++++++------------- 2 files changed, 254 insertions(+), 128 deletions(-) diff --git a/lang/m2/libm2/Mathlib.def b/lang/m2/libm2/Mathlib.def index ba595831c..5dcb0a10e 100644 --- a/lang/m2/libm2/Mathlib.def +++ b/lang/m2/libm2/Mathlib.def @@ -8,19 +8,27 @@ DEFINITION MODULE Mathlib; (* Some mathematical constants: *) CONST - (* From: Handbook of Mathematical Functions - Edited by M. Abramowitz and I.A. Stegun - National Bureau of Standards - Applied Mathematics Series 55 + (* From: Computer Approximations + Hart, Cheney, e.a. + The SIAM Series in Applied Mathematics + John Wiley & Sons, INC. New York London Sydney, 1968 *) - pi = 3.141592653589793238462643; - twicepi = 6.283185307179586476925286; - halfpi = 1.570796326794896619231322; - quartpi = 0.785398163397448309615661; - e = 2.718281828459045235360287; - ln2 = 0.693147180559945309417232; - ln10 = 2.302585092994045684017992; + pi = 3.14159265358979323846264338327950288; + twicepi = 6.28318530717958647692528676655900576; + halfpi = 1.57079632679489661923132169163975144; + quartpi = 0.78539816339744830961566084581987572; + e = 2.71828182845904523536028747135266250; + ln2 = 0.69314718055994530941723212145817657; + 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 *) diff --git a/lang/m2/libm2/Mathlib.mod b/lang/m2/libm2/Mathlib.mod index 7b84f89d6..86fc58fea 100644 --- a/lang/m2/libm2/Mathlib.mod +++ b/lang/m2/libm2/Mathlib.mod @@ -14,15 +14,11 @@ IMPLEMENTATION MODULE Mathlib; FROM EM IMPORT FIF, FEF; 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 OneRadianInDegrees = 57.295779513082320876798155D; OneDegreeInRadians = 0.017453292519943295769237D; + Sqrt2 = 1.41421356237309504880168872420969808D; + OneOverSqrt2 = 0.70710678118654752440084436210484904D; (* basic functions *) @@ -32,8 +28,7 @@ IMPLEMENTATION MODULE Mathlib; END pow; PROCEDURE longpow(x: LONGREAL; i: INTEGER): LONGREAL; - VAR - val: LONGREAL; + VAR val: LONGREAL; ri: LONGREAL; BEGIN ri := FLOATD(i); @@ -93,7 +88,7 @@ IMPLEMENTATION MODULE Mathlib; temp := temp / 2.0D; exp := exp + 2; END; - FOR i := 0 TO 4 DO + FOR i := 0 TO 5 DO temp := 0.5D*(temp + x/temp); END; RETURN temp; @@ -105,42 +100,50 @@ IMPLEMENTATION MODULE Mathlib; END exp; PROCEDURE longexp(x: LONGREAL): LONGREAL; - (* - * n = floor(x / ln2), d = x / ln2 - n - * exp(x) = exp((x / ln2) * ln2) = exp((n + d) * ln2) = - * exp(n * ln2) * exp(d * ln2) = 2 ** n * exp(d * ln2) - *) + (* 2**x = (Q(x*x)+x*P(x*x))/(Q(x*x)-x*P(x*x)) for x in [0,0.5] *) + (* Hart & Cheney #1069 *) CONST - a1 = -0.9999999995D; - a2 = 0.4999999206D; - a3 = -0.1666653019D; - a4 = 0.0416573475D; - a5 = -0.0083013598D; - a6 = 0.0013298820D; - a7 = -0.0001413161D; + p0 = 0.2080384346694663001443843411D+07; + p1 = 0.3028697169744036299076048876D+05; + p2 = 0.6061485330061080841615584556D+02; + q0 = 0.6002720360238832528230907598D+07; + q1 = 0.3277251518082914423057964422D+06; + q2 = 0.1749287689093076403844945335D+04; + q3 = 0.1000000000000000000000000000D+01; + VAR neg: BOOLEAN; - polval: LONGREAL; + xPxx, Qxx: LONGREAL; n: LONGREAL; n1 : INTEGER; + xsq : LONGREAL; + large: BOOLEAN; BEGIN neg := x < 0.0D; IF neg THEN x := -x; END; - x := FIF(x, 1.0D/LONG(ln2), n) * LONG(ln2); - polval := 1.0D /(1.0D + x*(a1+x*(a2+x*(a3+x*(a4+x*(a5+x*(a6+x*a7))))))); + x := FIF(x/longln2, 1.0D, n); + 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); WHILE n1 >= 16 DO - polval := polval * 65536.0D; + x := x * 65536.0D; n1 := n1 - 16; END; WHILE n1 > 0 DO - polval := polval * 2.0D; + x := x * 2.0D; DEC(n1); END; - IF neg THEN RETURN 1.0D/polval; END; - RETURN polval; + IF neg THEN RETURN 1.0D/x; END; + RETURN x; END longexp; PROCEDURE ln(x: REAL): REAL; (* natural log *) @@ -149,18 +152,23 @@ IMPLEMENTATION MODULE Mathlib; END ln; 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 - a1 = 0.9999964239D; - a2 = -0.4998741238D; - a3 = 0.3317990258D; - a4 = -0.2407338084D; - a5 = 0.1676540711D; - a6 = -0.0953293897D; - a7 = 0.0360884937D; - a8 = -0.0064535442D; + p0 = 0.7504094990777122217455611007D+02; + p1 = -0.1345669115050430235318253537D+03; + p2 = 0.7413719213248602512779336470D+02; + p3 = -0.1277249755012330819984385000D+02; + p4 = 0.3327108381087686938144000000D+00; + q0 = 0.3752047495388561108727775374D+02; + q1 = -0.7979028073715004879439951583D+02; + q2 = 0.5616126132118257292058560360D+02; + q3 = -0.1450868091858082685362325000D+02; + q4 = 0.1000000000000000000000000000D+01; VAR exp: INTEGER; - polval: LONGREAL; + z, zsq: LONGREAL; BEGIN IF x <= 0.0D THEN @@ -168,13 +176,15 @@ IMPLEMENTATION MODULE Mathlib; HALT END; x := FEF(x, exp); - WHILE x < 1.0D DO + WHILE x < OneOverSqrt2 DO x := x + x; DEC(exp); END; - x := x - 1.0D; - polval := x*(a1+x*(a2+x*(a3+x*(a4+x*(a5+x*(a6+x*(a7+a8*x))))))); - RETURN polval + FLOATD(exp) * LONG(ln2); + z := (x - 1.0D) / (x + 1.0D); + zsq := z*z; + 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; 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 *) BEGIN - RETURN longln(x)/LONG(ln10); + RETURN longln(x)/longln10; END longlog; (* trigonometric functions; arguments in radians *) @@ -194,34 +204,80 @@ IMPLEMENTATION MODULE Mathlib; RETURN SHORT(longsin(LONG(x))); 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 - a2 = -0.1666666664D; - a4 = 0.0083333315D; - a6 = -0.0001984090D; - a8 = 0.0000027526D; - a10 = -0.0000000239D; + p0 = 0.4857791909822798473837058825D+10; + p1 = -0.1808816670894030772075877725D+10; + p2 = 0.1724314784722489597789244188D+09; + p3 = -0.6351331748520454245913645971D+07; + 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 - xsqr: LONGREAL; - neg: BOOLEAN; + xsq, x1, x2, n : LONGREAL; + t : INTEGER; BEGIN - neg := FALSE; IF x < 0.0D THEN - neg := TRUE; + INC(quadrant, 2); x := -x; END; - x := FIF(x, 1.0D / LONG(twicepi), (* dummy *) xsqr) * LONG(twicepi); - IF x >= LONG(pi) THEN - neg := NOT neg; - x := x - LONG(pi); + IF longhalfpi - x = longhalfpi THEN + CASE quadrant OF + | 0,2: + RETURN 0.0D; + | 1: + RETURN 1.0D; + | 3: + RETURN -1.0D; + END; END; - IF x > LONG(halfpi) THEN - x := LONG(pi) - x; + IF x >= longtwicepi THEN + 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; - xsqr := x * x; - x := x * (1.0D + xsqr*(a2+xsqr*(a4+xsqr*(a6+xsqr*(a8+xsqr*a10))))); - IF neg THEN RETURN -x; END; - RETURN x; + x := x / longhalfpi; + t := TRUNC(x); + x := x - FLOATD(t); + 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; PROCEDURE cos(x: REAL): REAL; @@ -230,30 +286,9 @@ IMPLEMENTATION MODULE Mathlib; END cos; 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 - neg := FALSE; IF x < 0.0D THEN x := -x; END; - x := FIF(x, 1.0D / LONG(twicepi), (* dummy *) xsqr) * LONG(twicepi); - 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; + RETURN sinus(x, 1); END longcos; PROCEDURE tan(x: REAL): REAL; @@ -277,24 +312,31 @@ IMPLEMENTATION MODULE Mathlib; RETURN SHORT(longarcsin(LONG(x))); END arcsin; - PROCEDURE longarcsin(x: LONGREAL): LONGREAL; - CONST - a0 = 1.5707963050D; - a1 = -0.2145988016D; - a2 = 0.0889789874D; - a3 = -0.0501743046D; - a4 = 0.0308918810D; - a5 = -0.0170881256D; - a6 = 0.0066700901D; - a7 = -0.0012624911D; + PROCEDURE arcsincos(x: LONGREAL; cosfl: BOOLEAN): LONGREAL; + VAR + negative : BOOLEAN; BEGIN - IF x < 0.0D THEN x := -x; END; + negative := x <= 0.0D; + IF negative THEN x := -x; END; IF x > 1.0D THEN - Message("arcsin: argument > 1"); + Message("arcsin or arccos: argument > 1"); HALT END; - RETURN LONG(halfpi) - - longsqrt(1.0D - x)*(a0+x*(a1+x*(a2+x*(a3+x*(a4+x*(a5+x*(a6+x*a7))))))); + IF x = 1.0D THEN + 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; PROCEDURE arccos(x: REAL): REAL; @@ -304,7 +346,7 @@ IMPLEMENTATION MODULE Mathlib; PROCEDURE longarccos(x: LONGREAL): LONGREAL; BEGIN - RETURN LONG(halfpi) - longarcsin(x); + RETURN arcsincos(x, TRUE); END longarccos; PROCEDURE arctan(x: REAL): REAL; @@ -312,35 +354,109 @@ IMPLEMENTATION MODULE Mathlib; RETURN SHORT(longarctan(LONG(x))); 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; + (* 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 - a2 = -0.3333314528D; - a4 = 0.1999355085D; - a6 = -0.1420889944D; - a8 = 0.1065626393D; - a10 = -0.0752896400D; - a12 = 0.0429096318D; - a14 = -0.0161657367D; - a16 = 0.0028662257D; + p0 = 0.7698297257888171026986294745D+03; + p1 = 0.1557282793158363491416585283D+04; + p2 = 0.1033384651675161628243434662D+04; + p3 = 0.2485841954911840502660889866D+03; + p4 = 0.1566564964979791769948970100D+02; + q0 = 0.7698297257888171026986294911D+03; + q1 = 0.1813892701754635858982709369D+04; + q2 = 0.1484049607102276827437401170D+04; + q3 = 0.4904645326203706217748848797D+03; + q4 = 0.5593479839280348664778328000D+02; + q5 = 0.1000000000000000000000000000D+01; VAR xsqr: LONGREAL; - rev: BOOLEAN; neg: BOOLEAN; - BEGIN - rev := FALSE; + i: INTEGER; + BEGIN + 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; IF x < 0.0D THEN neg := TRUE; x := -x; END; - IF x > 1.0D THEN - rev := TRUE; - x := 1.0D / x; + i := 0; + WHILE (i <= 3) AND (x <= precomp[i].X) DO + INC(i); END; - xsqr := x * x; - x := x * (1.0D + xsqr*(a2+xsqr*(a4+xsqr*(a6+xsqr*(a8+xsqr*(a10+xsqr*(a12+xsqr*(a14+xsqr*a16)))))))); - IF rev THEN - x := LONG(quartpi) - x; + IF (i # 0) THEN + WITH precomp[i] DO + x := arctan + longarctan(OneOverXn-OneOverXnSquarePlusone/(OneOverXn+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; IF neg THEN RETURN -x; END; RETURN x; @@ -452,4 +568,6 @@ IMPLEMENTATION MODULE Mathlib; RETURN x * OneDegreeInRadians; END longDegreeToRadian; +BEGIN + arctaninit := FALSE; END Mathlib. -- 2.34.1