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 *)
END pow;
PROCEDURE longpow(x: LONGREAL; i: INTEGER): LONGREAL;
- VAR
- val: LONGREAL;
+ VAR val: LONGREAL;
ri: LONGREAL;
BEGIN
ri := FLOATD(i);
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;
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 *)
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
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 *)
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 *)
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;
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;
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;
PROCEDURE longarccos(x: LONGREAL): LONGREAL;
BEGIN
- RETURN LONG(halfpi) - longarcsin(x);
+ RETURN arcsincos(x, TRUE);
END longarccos;
PROCEDURE arctan(x: REAL): REAL;
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;
RETURN x * OneDegreeInRadians;
END longDegreeToRadian;
+BEGIN
+ arctaninit := FALSE;
END Mathlib.
public / ack.git / commitdiff