123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239240241242243244245246247248249250251252253254255256257258259260261262263264265266267268269270271272273274275276277278279280281282283284285286287288289290291292293294295296297298299300301302303304305306307308309310311312313314315316317318319320321322323324325326327328329330331332333334335336337338339340341342343344345346347348349350351352353354355356357358359360361362363364365366367368369370371372373374375376377378379380381382383384385386387388389390391392393394395396397398399400401402403404405406407408409410411412413414415416417418419420421422423424425426427428429430431432433434435436437438439440441442443444445446447448449450451452453454455456457458459460461462463464465466467468469470471472473474475476477478479480481482483484485486487488489490491492493494495496497498499500501502503504505506507508509510511512513514515516517518519520521522523524525526527528529530531532533534535536537538539540541542543544545546547548549550551552553554555556557558559560561562563564565566567568569570571572573574575576 |
- (*
- (c) copyright 1988 by the Vrije Universiteit, Amsterdam, The Netherlands.
- See the copyright notice in the ACK home directory, in the file "Copyright".
- *)
- (*$R-*)
- IMPLEMENTATION MODULE Mathlib;
- (*
- Module: Mathematical functions
- Author: Ceriel J.H. Jacobs
- Version: $Id$
- *)
- FROM EM IMPORT FIF, FEF;
- FROM Traps IMPORT Message;
- CONST
- OneRadianInDegrees = 57.295779513082320876798155D;
- OneDegreeInRadians = 0.017453292519943295769237D;
- OneOverSqrt2 = 0.70710678118654752440084436210484904D;
- (* basic functions *)
- PROCEDURE pow(x: REAL; i: INTEGER): REAL;
- BEGIN
- RETURN SHORT(longpow(LONG(x), i));
- END pow;
- PROCEDURE longpow(x: LONGREAL; i: INTEGER): LONGREAL;
- VAR val: LONGREAL;
- ri: LONGREAL;
- BEGIN
- ri := FLOATD(i);
- IF x < 0.0D THEN
- val := longexp(longln(-x) * ri);
- IF ODD(i) THEN RETURN -val;
- ELSE RETURN val;
- END;
- ELSIF x = 0.0D THEN
- RETURN 0.0D;
- ELSE
- RETURN longexp(longln(x) * ri);
- END;
- END longpow;
- PROCEDURE sqrt(x: REAL): REAL;
- BEGIN
- RETURN SHORT(longsqrt(LONG(x)));
- END sqrt;
- PROCEDURE longsqrt(x: LONGREAL): LONGREAL;
- VAR
- temp: LONGREAL;
- exp, i: INTEGER;
- BEGIN
- IF x <= 0.0D THEN
- IF x < 0.0D THEN
- Message("sqrt: negative argument");
- HALT
- END;
- RETURN 0.0D;
- END;
- temp := FEF(x,exp);
- (*
- * NOTE
- * this wont work on 1's comp
- *)
- IF ODD(exp) THEN
- temp := 2.0D * temp;
- DEC(exp);
- END;
- temp := 0.5D*(1.0D + temp);
- WHILE exp > 28 DO
- temp := temp * 16384.0D;
- exp := exp - 28;
- END;
- WHILE exp < -28 DO
- temp := temp / 16384.0D;
- exp := exp + 28;
- END;
- WHILE exp >= 2 DO
- temp := temp * 2.0D;
- exp := exp - 2;
- END;
- WHILE exp <= -2 DO
- temp := temp / 2.0D;
- exp := exp + 2;
- END;
- FOR i := 0 TO 5 DO
- temp := 0.5D*(temp + x/temp);
- END;
- RETURN temp;
- 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;
- BEGIN
- RETURN SHORT(longexp(LONG(x)));
- END exp;
- PROCEDURE longexp(x: LONGREAL): LONGREAL;
- (* Algorithm and coefficients from:
- "Software manual for the elementary functions"
- by W.J. Cody and W. Waite, Prentice-Hall, 1980
- *)
- CONST
- p0 = 0.25000000000000000000D+00;
- p1 = 0.75753180159422776666D-02;
- p2 = 0.31555192765684646356D-04;
- q0 = 0.50000000000000000000D+00;
- q1 = 0.56817302698551221787D-01;
- q2 = 0.63121894374398503557D-03;
- q3 = 0.75104028399870046114D-06;
- VAR
- neg: BOOLEAN;
- n: INTEGER;
- xn, g, x1, x2: LONGREAL;
- BEGIN
- neg := x < 0.0D;
- IF neg THEN
- x := -x;
- END;
- n := TRUNC(x/longln2 + 0.5D);
- xn := FLOATD(n);
- x1 := FLOATD(TRUNCD(x));
- x2 := x - x1;
- g := ((x1 - xn * 0.693359375D)+x2) - xn * (-2.1219444005469058277D-4);
- IF neg THEN
- g := -g;
- n := -n;
- END;
- xn := g*g;
- x := g*((p2*xn+p1)*xn+p0);
- INC(n);
- RETURN ldexp(0.5D + x/((((q3*xn+q2)*xn+q1)*xn+q0) - x), n);
- END longexp;
- PROCEDURE ln(x: REAL): REAL; (* natural log *)
- BEGIN
- RETURN SHORT(longln(LONG(x)));
- END ln;
- PROCEDURE longln(x: LONGREAL): LONGREAL; (* natural log *)
- (* Algorithm and coefficients from:
- "Software manual for the elementary functions"
- by W.J. Cody and W. Waite, Prentice-Hall, 1980
- *)
- CONST
- p0 = -0.64124943423745581147D+02;
- p1 = 0.16383943563021534222D+02;
- p2 = -0.78956112887491257267D+00;
- q0 = -0.76949932108494879777D+03;
- q1 = 0.31203222091924532844D+03;
- q2 = -0.35667977739034646171D+02;
- q3 = 1.0D;
- VAR
- exp: INTEGER;
- z, znum, zden, w: LONGREAL;
- BEGIN
- IF x <= 0.0D THEN
- Message("ln: argument <= 0");
- HALT
- END;
- x := FEF(x, exp);
- IF x > OneOverSqrt2 THEN
- 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);
- END;
- z := znum / zden;
- w := z * z;
- x := z + z * w * (((p2*w+p1)*w+p0)/(((q3*w+q2)*w+q1)*w+q0));
- z := FLOATD(exp);
- x := x + z * (-2.121944400546905827679D-4);
- RETURN x + z * 0.693359375D;
- END longln;
- PROCEDURE log(x: REAL): REAL; (* log with base 10 *)
- BEGIN
- RETURN SHORT(longlog(LONG(x)));
- END log;
- PROCEDURE longlog(x: LONGREAL): LONGREAL; (* log with base 10 *)
- BEGIN
- RETURN longln(x)/longln10;
- END longlog;
- (* trigonometric functions; arguments in radians *)
- PROCEDURE sin(x: REAL): REAL;
- BEGIN
- RETURN SHORT(longsin(LONG(x)));
- END sin;
- PROCEDURE sinus(x: LONGREAL; cosflag: BOOLEAN) : LONGREAL;
- (* Algorithm and coefficients from:
- "Software manual for the elementary functions"
- by W.J. Cody and W. Waite, Prentice-Hall, 1980
- *)
- CONST
- r0 = -0.16666666666666665052D+00;
- r1 = 0.83333333333331650314D-02;
- r2 = -0.19841269841201840457D-03;
- r3 = 0.27557319210152756119D-05;
- r4 = -0.25052106798274584544D-07;
- r5 = 0.16058936490371589114D-09;
- r6 = -0.76429178068910467734D-12;
- r7 = 0.27204790957888846175D-14;
- A1 = 3.1416015625D;
- A2 = -8.908910206761537356617D-6;
- VAR
- x1, x2, y : LONGREAL;
- neg : BOOLEAN;
- BEGIN
- IF x < 0.0D THEN
- neg := TRUE;
- x := -x
- ELSE neg := FALSE
- END;
- IF cosflag THEN
- neg := FALSE;
- y := longhalfpi + x
- ELSE
- y := x
- END;
- y := y / longpi + 0.5D;
- IF FIF(y, 1.0D, y) < 0.0D THEN ; END;
- IF FIF(y, 0.5D, x1) # 0.0D THEN neg := NOT neg END;
- IF cosflag THEN y := y - 0.5D END;
- x2 := FIF(x, 1.0, x1);
- x := x1 - y * A1;
- x := x + x2;
- x := x - y * A2;
- IF x < 0.0D THEN
- neg := NOT neg;
- x := -x
- END;
- y := x * x;
- x := x + x * y * (((((((r7*y+r6)*y+r5)*y+r4)*y+r3)*y+r2)*y+r1)*y+r0);
- IF neg THEN RETURN -x END;
- RETURN x;
- END sinus;
- PROCEDURE longsin(x: LONGREAL): LONGREAL;
- BEGIN
- RETURN sinus(x, FALSE);
- END longsin;
- PROCEDURE cos(x: REAL): REAL;
- BEGIN
- RETURN SHORT(longcos(LONG(x)));
- END cos;
- PROCEDURE longcos(x: LONGREAL): LONGREAL;
- BEGIN
- IF x < 0.0D THEN x := -x; END;
- RETURN sinus(x, TRUE);
- END longcos;
- PROCEDURE tan(x: REAL): REAL;
- BEGIN
- RETURN SHORT(longtan(LONG(x)));
- END tan;
- PROCEDURE longtan(x: LONGREAL): 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
- negative := x < 0.0D;
- y := x / longhalfpi + 0.5D;
- (* Use extended precision to calculate reduced argument.
- Here we used 12 bits of the mantissa for a1.
- 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;
- PROCEDURE arcsin(x: REAL): REAL;
- BEGIN
- RETURN SHORT(longarcsin(LONG(x)));
- END arcsin;
- PROCEDURE arcsincos(x: LONGREAL; cosfl: BOOLEAN): LONGREAL;
- 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;
- big: BOOLEAN;
- g: LONGREAL;
- BEGIN
- negative := x < 0.0D;
- IF negative THEN x := -x; END;
- IF x > 0.5D THEN
- big := TRUE;
- IF x > 1.0D THEN
- Message("arcsin or arccos: argument > 1");
- HALT
- END;
- g := 0.5D - 0.5D * x;
- x := -longsqrt(g);
- x := x + x;
- ELSE
- big := FALSE;
- g := x * x;
- END;
- x := x + x * g *
- ((((p4*g+p3)*g+p2)*g+p1)*g+p0)/(((((q5*g+q4)*g+q3)*g+q2)*g+q1)*g+q0);
- 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;
- IF negative AND NOT cosfl THEN x := -x END;
- RETURN x;
- END arcsincos;
- PROCEDURE longarcsin(x: LONGREAL): LONGREAL;
- BEGIN
- RETURN arcsincos(x, FALSE);
- END longarcsin;
- PROCEDURE arccos(x: REAL): REAL;
- BEGIN
- RETURN SHORT(longarccos(LONG(x)));
- END arccos;
- PROCEDURE longarccos(x: LONGREAL): LONGREAL;
- BEGIN
- RETURN arcsincos(x, TRUE);
- END longarccos;
- PROCEDURE arctan(x: REAL): REAL;
- BEGIN
- RETURN SHORT(longarctan(LONG(x)));
- END arctan;
- VAR A: ARRAY[0..3] OF LONGREAL;
- arctaninit: BOOLEAN;
- PROCEDURE longarctan(x: LONGREAL): LONGREAL;
- (* Algorithm and coefficients from:
- "Software manual for the elementary functions"
- by W.J. Cody and W. Waite, Prentice-Hall, 1980
- *)
- CONST
- p0 = -0.13688768894191926929D+02;
- p1 = -0.20505855195861651981D+02;
- p2 = -0.84946240351320683534D+01;
- p3 = -0.83758299368150059274D+00;
- q0 = 0.41066306682575781263D+02;
- q1 = 0.86157349597130242515D+02;
- q2 = 0.59578436142597344465D+02;
- q3 = 0.15024001160028576121D+02;
- q4 = 1.0D;
- VAR
- g: LONGREAL;
- neg: BOOLEAN;
- n: INTEGER;
- BEGIN
- IF NOT arctaninit THEN
- arctaninit := TRUE;
- A[0] := 0.0D;
- A[1] := 0.52359877559829887307710723554658381D; (* p1/6 *)
- A[2] := longhalfpi;
- A[3] := 1.04719755119659774615421446109316763D; (* pi/3 *)
- END;
- neg := FALSE;
- IF x < 0.0D THEN
- neg := TRUE;
- x := -x;
- END;
- IF x > 1.0D THEN
- x := 1.0D/x;
- n := 2
- ELSE
- n := 0
- 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;
- RETURN x;
- END longarctan;
- (* 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;
- BEGIN
- RETURN SHORT(longsinh(LONG(x)));
- END sinh;
- PROCEDURE longsinh(x: LONGREAL): LONGREAL;
- VAR expx: LONGREAL;
- BEGIN
- expx := longexp(x);
- RETURN (expx - 1.0D/expx)/2.0D;
- END longsinh;
- PROCEDURE cosh(x: REAL): REAL;
- BEGIN
- RETURN SHORT(longcosh(LONG(x)));
- END cosh;
- PROCEDURE longcosh(x: LONGREAL): LONGREAL;
- VAR expx: LONGREAL;
- BEGIN
- expx := longexp(x);
- RETURN (expx + 1.0D/expx)/2.0D;
- END longcosh;
- PROCEDURE tanh(x: REAL): REAL;
- BEGIN
- RETURN SHORT(longtanh(LONG(x)));
- END tanh;
- PROCEDURE longtanh(x: LONGREAL): LONGREAL;
- VAR expx: LONGREAL;
- BEGIN
- expx := longexp(x);
- RETURN (expx - 1.0D/expx) / (expx + 1.0D/expx);
- END longtanh;
- PROCEDURE arcsinh(x: REAL): REAL;
- BEGIN
- RETURN SHORT(longarcsinh(LONG(x)));
- END arcsinh;
- PROCEDURE longarcsinh(x: LONGREAL): LONGREAL;
- VAR neg: BOOLEAN;
- BEGIN
- neg := FALSE;
- IF x < 0.0D THEN
- neg := TRUE;
- x := -x;
- END;
- x := longln(x + longsqrt(x*x+1.0D));
- IF neg THEN RETURN -x; END;
- RETURN x;
- END longarcsinh;
- PROCEDURE arccosh(x: REAL): REAL;
- BEGIN
- RETURN SHORT(longarccosh(LONG(x)));
- END arccosh;
- PROCEDURE longarccosh(x: LONGREAL): LONGREAL;
- BEGIN
- IF x < 1.0D THEN
- Message("arccosh: argument < 1");
- HALT
- END;
- RETURN longln(x + longsqrt(x*x - 1.0D));
- END longarccosh;
- PROCEDURE arctanh(x: REAL): REAL;
- BEGIN
- RETURN SHORT(longarctanh(LONG(x)));
- END arctanh;
- PROCEDURE longarctanh(x: LONGREAL): LONGREAL;
- BEGIN
- IF (x <= -1.0D) OR (x >= 1.0D) THEN
- Message("arctanh: ABS(argument) >= 1");
- HALT
- END;
- RETURN longln((1.0D + x)/(1.0D - x)) / 2.0D;
- END longarctanh;
- (* conversions *)
- PROCEDURE RadianToDegree(x: REAL): REAL;
- BEGIN
- RETURN SHORT(longRadianToDegree(LONG(x)));
- END RadianToDegree;
- PROCEDURE longRadianToDegree(x: LONGREAL): LONGREAL;
- BEGIN
- RETURN x * OneRadianInDegrees;
- END longRadianToDegree;
- PROCEDURE DegreeToRadian(x: REAL): REAL;
- BEGIN
- RETURN SHORT(longDegreeToRadian(LONG(x)));
- END DegreeToRadian;
- PROCEDURE longDegreeToRadian(x: LONGREAL): LONGREAL;
- BEGIN
- RETURN x * OneDegreeInRadians;
- END longDegreeToRadian;
- BEGIN
- arctaninit := FALSE;
- END Mathlib.
|