From 86c5c56a387c72efc19b1b8cb1e9eefd046bed97 Mon Sep 17 00:00:00 2001 From: ceriel Date: Wed, 27 May 1987 10:05:01 +0000 Subject: [PATCH] Added Mathlib; MathLib0 now uses Mathlib --- lang/m2/libm2/.distr | 1 + lang/m2/libm2/FIFFEF.def | 5 +- lang/m2/libm2/FIFFEF.e | 45 ++-- lang/m2/libm2/LIST | 5 +- lang/m2/libm2/Makefile | 3 +- lang/m2/libm2/MathLib0.def | 4 + lang/m2/libm2/MathLib0.mod | 305 +----------------------- lang/m2/libm2/Mathlib.def | 66 ++++++ lang/m2/libm2/Mathlib.mod | 443 +++++++++++++++++++++++++++++++++++ lang/m2/libm2/RealConver.mod | 65 +++-- 10 files changed, 591 insertions(+), 351 deletions(-) create mode 100644 lang/m2/libm2/Mathlib.def create mode 100644 lang/m2/libm2/Mathlib.mod diff --git a/lang/m2/libm2/.distr b/lang/m2/libm2/.distr index a91543de3..5c9d5afa5 100644 --- a/lang/m2/libm2/.distr +++ b/lang/m2/libm2/.distr @@ -6,6 +6,7 @@ Conversion.def FIFFEF.def InOut.def Makefile +Mathlib.def MathLib0.def Processes.def RealInOut.def diff --git a/lang/m2/libm2/FIFFEF.def b/lang/m2/libm2/FIFFEF.def index 922332d12..ea49df35a 100644 --- a/lang/m2/libm2/FIFFEF.def +++ b/lang/m2/libm2/FIFFEF.def @@ -1,11 +1,12 @@ +(*$Foreign*) DEFINITION MODULE FIFFEF; - PROCEDURE FIF(arg1, arg2: REAL; VAR intres: REAL) : REAL; + PROCEDURE FIF(arg1, arg2: LONGREAL; VAR intres: LONGREAL) : LONGREAL; (* multiplies arg1 and arg2, and returns the integer part of the result in "intres" and the fraction part as the function result. *) - PROCEDURE FEF(arg: REAL; VAR exp: INTEGER) : REAL; + PROCEDURE FEF(arg: LONGREAL; VAR exp: INTEGER) : LONGREAL; (* splits "arg" in mantissa and a base-2 exponent. The mantissa is returned, and the exponent is left in "exp". *) diff --git a/lang/m2/libm2/FIFFEF.e b/lang/m2/libm2/FIFFEF.e index b7588e6a9..849cfc2b7 100644 --- a/lang/m2/libm2/FIFFEF.e +++ b/lang/m2/libm2/FIFFEF.e @@ -2,50 +2,45 @@ mes 2,EM_WSIZE,EM_PSIZE #define ARG1 0 -#define ARG2 EM_FSIZE -#define IRES 2*EM_FSIZE +#define ARG2 EM_DSIZE +#define IRES 2*EM_DSIZE -; FIFFEF_FIF is called with three parameters: +; FIF is called with three parameters: ; - address of integer part result (IRES) ; - float two (ARG2) ; - float one (ARG1) -; and returns an EM_FSIZE-byte floating point number +; and returns an EM_DSIZE-byte floating point number ; Definition: -; PROCEDURE FIF(ARG1, ARG2: REAL; VAR IRES: REAL) : REAL; +; PROCEDURE FIF(ARG1, ARG2: LONGREAL; VAR IRES: LONGREAL) : LONGREAL; - exp $FIFFEF_FIF - pro $FIFFEF_FIF,0 + exp $FIF + pro $FIF,0 lal 0 - loi 2*EM_FSIZE - fif EM_FSIZE + loi 2*EM_DSIZE + fif EM_DSIZE lal IRES loi EM_PSIZE - sti EM_FSIZE - ret EM_FSIZE + sti EM_DSIZE + ret EM_DSIZE end ? #define FARG 0 -#define ERES EM_FSIZE +#define ERES EM_DSIZE -; FIFFEF_FEF is called with two parameters: +; FEF is called with two parameters: ; - address of base 2 exponent result (ERES) ; - floating point number to be split (FARG) -; and returns an EM_FSIZE-byte floating point number (the mantissa) +; and returns an EM_DSIZE-byte floating point number (the mantissa) ; Definition: -; PROCEDURE FEF(FARG: REAL; VAR ERES: integer): REAL; +; PROCEDURE FEF(FARG: LONGREAL; VAR ERES: integer): LONGREAL; - exp $FIFFEF_FEF - pro $FIFFEF_FEF,0 + exp $FEF + pro $FEF,0 lal FARG - loi EM_FSIZE - fef EM_FSIZE + loi EM_DSIZE + fef EM_DSIZE lal ERES loi EM_PSIZE sti EM_WSIZE - ret EM_FSIZE - end ? - - exp $FIFFEF - pro $FIFFEF,0 - ret 0 + ret EM_DSIZE end ? diff --git a/lang/m2/libm2/LIST b/lang/m2/libm2/LIST index 05d106755..9e2d8a06c 100644 --- a/lang/m2/libm2/LIST +++ b/lang/m2/libm2/LIST @@ -1,18 +1,19 @@ tail_m2.a +RealInOut.mod InOut.mod Terminal.mod TTY.mod ASCII.mod -FIFFEF.e MathLib0.mod +Mathlib.mod Processes.mod RealConver.mod -RealInOut.mod Storage.mod Conversion.mod Semaphores.mod random.mod Strings.mod +FIFFEF.e Arguments.c catch.c hol0.e diff --git a/lang/m2/libm2/Makefile b/lang/m2/libm2/Makefile index 1b2dbf2f1..77b7bc94b 100644 --- a/lang/m2/libm2/Makefile +++ b/lang/m2/libm2/Makefile @@ -4,7 +4,8 @@ DEFDIR = $(HOME)/lib/m2 SOURCES = ASCII.def FIFFEF.def MathLib0.def Processes.def \ RealInOut.def Storage.def Arguments.def Conversion.def \ random.def Semaphores.def Unix.def RealConver.def \ - Strings.def InOut.def Terminal.def TTY.def + Strings.def InOut.def Terminal.def TTY.def \ + Mathlib.def all: diff --git a/lang/m2/libm2/MathLib0.def b/lang/m2/libm2/MathLib0.def index cbae3a382..985a6da44 100644 --- a/lang/m2/libm2/MathLib0.def +++ b/lang/m2/libm2/MathLib0.def @@ -1,4 +1,8 @@ DEFINITION MODULE MathLib0; +(* + Exists for compatibility. + A more elaborate math lib can be found in Mathlib.def +*) PROCEDURE sqrt(x : REAL) : REAL; diff --git a/lang/m2/libm2/MathLib0.mod b/lang/m2/libm2/MathLib0.mod index 93c8a7864..d48833c47 100644 --- a/lang/m2/libm2/MathLib0.mod +++ b/lang/m2/libm2/MathLib0.mod @@ -1,326 +1,35 @@ IMPLEMENTATION MODULE MathLib0; -(* Rewritten in Modula-2. - The originals came from the Pascal runtime library. -*) -FROM FIFFEF IMPORT FIF, FEF; - -CONST - HUGE = 1.701411733192644270E38; - -PROCEDURE sinus(arg: REAL; quad: INTEGER): REAL; - -(* Coefficients for sin/cos are #3370 from Hart & Cheney (18.80D). -*) -CONST - twoopi = 0.63661977236758134308; - p0 = 0.1357884097877375669092680E8; - p1 = -0.4942908100902844161158627E7; - p2 = 0.4401030535375266501944918E6; - p3 = -0.1384727249982452873054457E5; - p4 = 0.1459688406665768722226959E3; - q0 = 0.8644558652922534429915149E7; - q1 = 0.4081792252343299749395779E6; - q2 = 0.9463096101538208180571257E4; - q3 = 0.1326534908786136358911494E3; -VAR - e, f: REAL; - ysq: REAL; - x,y: REAL; - k: INTEGER; - temp1, temp2: REAL; -BEGIN - x := arg; - IF x < 0.0 THEN - x := -x; - quad := quad + 2; - END; - x := x*twoopi; (*underflow?*) - IF x>32764.0 THEN - y := FIF(x, 10.0, e); - e := e + FLOAT(quad); - temp1 := FIF(0.25, e, f); - quad := TRUNC(e - 4.0*f); - ELSE - k := TRUNC(x); - y := x - FLOAT(k); - quad := (quad + k) MOD 4; - END; - IF ODD(quad) THEN - y := 1.0-y; - END; - IF quad > 1 THEN - y := -y; - END; - - ysq := y*y; - temp1 := ((((p4*ysq+p3)*ysq+p2)*ysq+p1)*ysq+p0)*y; - temp2 := ((((ysq+q3)*ysq+q2)*ysq+q1)*ysq+q0); - RETURN temp1/temp2; -END sinus; + IMPORT Mathlib; PROCEDURE cos(arg: REAL): REAL; BEGIN - IF arg < 0.0 THEN - arg := -arg; - END; - RETURN sinus(arg, 1); + RETURN Mathlib.cos(arg); END cos; PROCEDURE sin(arg: REAL): REAL; BEGIN - RETURN sinus(arg, 0); + RETURN Mathlib.sin(arg); END sin; -(* - floating-point arctangent - - arctan returns the value of the arctangent of its - argument in the range [-pi/2,pi/2]. - - coefficients are #5077 from Hart & Cheney. (19.56D) -*) - -CONST - sq2p1 = 2.414213562373095048802E0; - sq2m1 = 0.414213562373095048802E0; - pio2 = 1.570796326794896619231E0; - pio4 = 0.785398163397448309615E0; - p4 = 0.161536412982230228262E2; - p3 = 0.26842548195503973794141E3; - p2 = 0.11530293515404850115428136E4; - p1 = 0.178040631643319697105464587E4; - p0 = 0.89678597403663861959987488E3; - q4 = 0.5895697050844462222791E2; - q3 = 0.536265374031215315104235E3; - q2 = 0.16667838148816337184521798E4; - q1 = 0.207933497444540981287275926E4; - q0 = 0.89678597403663861962481162E3; - -(* - xatan evaluates a series valid in the - range [-0.414...,+0.414...]. -*) - -PROCEDURE xatan(arg: REAL) : REAL; -VAR - argsq, value: REAL; -BEGIN - argsq := arg*arg; - value := ((((p4*argsq + p3)*argsq + p2)*argsq + p1)*argsq + p0); - value := value/(((((argsq + q4)*argsq + q3)*argsq + q2)*argsq + q1)*argsq + q0); - RETURN value*arg; -END xatan; - -PROCEDURE satan(arg: REAL): REAL; -BEGIN - IF arg < sq2m1 THEN - RETURN xatan(arg); - ELSIF arg > sq2p1 THEN - RETURN pio2 - xatan(1.0/arg); - ELSE - RETURN pio4 + xatan((arg-1.0)/(arg+1.0)); - END; -END satan; - -(* - atan makes its argument positive and - calls the inner routine satan. -*) - PROCEDURE arctan(arg: REAL): REAL; BEGIN - IF arg>0.0 THEN - RETURN satan(arg); - ELSE - RETURN -satan(-arg); - END; + RETURN Mathlib.arctan(arg); END arctan; -(* - sqrt returns the square root of its floating - point argument. Newton's method. -*) - PROCEDURE sqrt(arg: REAL): REAL; -VAR - x, temp: REAL; - exp, i: INTEGER; BEGIN - IF arg <= 0.0 THEN - IF arg < 0.0 THEN - (* ??? *) - ; - END; - RETURN 0.0; - END; - x := FEF(arg,exp); - (* - * NOTE - * this wont work on 1's comp - *) - IF ODD(exp) THEN - x := 2.0 * x; - DEC(exp); - END; - temp := 0.5*(1.0 + x); - - WHILE exp > 28 DO - temp := temp * 16384.0; - exp := exp - 28; - END; - WHILE exp < -28 DO - temp := temp / 16384.0; - exp := exp + 28; - END; - WHILE exp >= 2 DO - temp := temp * 2.0; - exp := exp - 2; - END; - WHILE exp <= -2 DO - temp := temp / 2.0; - exp := exp + 2; - END; - FOR i := 0 TO 4 DO - temp := 0.5*(temp + arg/temp); - END; - RETURN temp; + RETURN Mathlib.sqrt(arg); END sqrt; -(* - ln returns the natural logarithm of its floating - point argument. - - The coefficients are #2705 from Hart & Cheney. (19.38D) -*) PROCEDURE ln(arg: REAL): REAL; -CONST - log2 = 0.693147180559945309E0; - sqrto2 = 0.707106781186547524E0; - p0 = -0.240139179559210510E2; - p1 = 0.309572928215376501E2; - p2 = -0.963769093368686593E1; - p3 = 0.421087371217979714E0; - q0 = -0.120069589779605255E2; - q1 = 0.194809660700889731E2; - q2 = -0.891110902798312337E1; -VAR - x,z, zsq, temp: REAL; - exp: INTEGER; BEGIN - IF arg <= 0.0 THEN - (* ??? *) - RETURN -HUGE; - END; - x := FEF(arg,exp); - IF x 127 THEN - (* Too large. ??? *) - RETURN FLOAT(neg) * HUGE; - END; - IF exp < -127 THEN - RETURN 0.0; - END; - WHILE exp > 14 DO - fr := fr * 16384.0; - exp := exp - 14; - END; - WHILE exp < -14 DO - fr := fr / 16384.0; - exp := exp + 14; - END; - WHILE exp > 0 DO - fr := fr + fr; - DEC(exp); - END; - WHILE exp < 0 DO - fr := fr / 2.0; - INC(exp); - END; - RETURN FLOAT(neg) * fr; -END ldexp; - PROCEDURE exp(arg: REAL): REAL; -CONST - p0 = 0.2080384346694663001443843411E7; - p1 = 0.3028697169744036299076048876E5; - p2 = 0.6061485330061080841615584556E2; - q0 = 0.6002720360238832528230907598E7; - q1 = 0.3277251518082914423057964422E6; - q2 = 0.1749287689093076403844945335E4; - log2e = 1.4426950408889634073599247; - sqrt2 = 1.4142135623730950488016887; - maxf = 10000.0; -VAR - fract: REAL; - temp1, temp2, xsq: REAL; - ent: INTEGER; BEGIN - IF arg = 0.0 THEN - RETURN 1.0; - END; - IF arg < -maxf THEN - RETURN 0.0; - END; - IF arg > maxf THEN - (* result too large ??? *) - RETURN HUGE; - END; - arg := arg * log2e; - ent := TRUNC(floor(arg)); - fract := (arg-FLOAT(ent)) - 0.5; - xsq := fract*fract; - temp1 := ((p2*xsq+p1)*xsq+p0)*fract; - temp2 := ((xsq+q2)*xsq+q1)*xsq + q0; - RETURN ldexp(sqrt2*(temp2+temp1)/(temp2-temp1), ent); + RETURN Mathlib.exp(arg); END exp; PROCEDURE entier(x: REAL): INTEGER; diff --git a/lang/m2/libm2/Mathlib.def b/lang/m2/libm2/Mathlib.def new file mode 100644 index 000000000..979d1b696 --- /dev/null +++ b/lang/m2/libm2/Mathlib.def @@ -0,0 +1,66 @@ +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 + *) + + pi = 3.141592653589793238462643; + twicepi = 6.283185307179586476925286; + halfpi = 1.570796326794896619231322; + quartpi = 0.785398163397448309615661; + e = 2.718281828459045235360287; + ln2 = 0.693147180559945309417232; + ln10 = 2.302585092994045684017992; + + (* basic functions *) + + PROCEDURE pow(x: REAL; i: INTEGER): REAL; + + PROCEDURE sqrt(x: REAL): REAL; + + PROCEDURE exp(x: REAL): REAL; + + PROCEDURE ln(x: REAL): REAL; (* natural log *) + + PROCEDURE log(x: REAL): REAL; (* log with base 10 *) + + (* trigonometric functions; arguments in radians *) + + PROCEDURE sin(x: REAL): REAL; + + PROCEDURE cos(x: REAL): REAL; + + PROCEDURE tan(x: REAL): REAL; + + PROCEDURE arcsin(x: REAL): REAL; + + PROCEDURE arccos(x: REAL): REAL; + + PROCEDURE arctan(x: REAL): REAL; + + (* hyperbolic functions *) + + PROCEDURE sinh(x: REAL): REAL; + + PROCEDURE cosh(x: REAL): REAL; + + PROCEDURE tanh(x: REAL): REAL; + + PROCEDURE arcsinh(x: REAL): REAL; + + PROCEDURE arccosh(x: REAL): REAL; + + PROCEDURE arctanh(x: REAL): REAL; + + (* conversions *) + + PROCEDURE RadianToDegree(x: REAL): REAL; + + PROCEDURE DegreeToRadian(x: REAL): REAL; + +END Mathlib. diff --git a/lang/m2/libm2/Mathlib.mod b/lang/m2/libm2/Mathlib.mod new file mode 100644 index 000000000..e1268bb32 --- /dev/null +++ b/lang/m2/libm2/Mathlib.mod @@ -0,0 +1,443 @@ +IMPLEMENTATION MODULE Mathlib; + + FROM FIFFEF IMPORT FIF, FEF; + + (* 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; + + (* 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 + (* ??? *) + ; + 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 4 DO + temp := 0.5D*(temp + x/temp); + END; + RETURN temp; + END longsqrt; + + PROCEDURE exp(x: REAL): REAL; + BEGIN + RETURN SHORT(longexp(LONG(x))); + 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) + *) + CONST + a1 = -0.9999999995D; + a2 = 0.4999999206D; + a3 = -0.1666653019D; + a4 = 0.0416573475D; + a5 = -0.0083013598D; + a6 = 0.0013298820D; + a7 = -0.0001413161D; + VAR + neg: BOOLEAN; + polval: LONGREAL; + n: LONGREAL; + n1 : INTEGER; + 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))))))); + n1 := TRUNCD(n + 0.5D); + WHILE n1 >= 16 DO + polval := polval * 65536.0D; + n1 := n1 - 16; + END; + WHILE n1 > 0 DO + polval := polval * 2.0D; + DEC(n1); + END; + IF neg THEN RETURN 1.0D/polval; END; + RETURN polval; + END longexp; + + PROCEDURE ln(x: REAL): REAL; (* natural log *) + BEGIN + RETURN SHORT(longln(LONG(x))); + END ln; + + PROCEDURE longln(x: LONGREAL): LONGREAL; (* natural log *) + 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; + VAR + exp: INTEGER; + polval: LONGREAL; + + BEGIN + IF x <= 0.0D THEN + (* ??? *) + RETURN 0.0D; + END; + x := FEF(x, exp); + WHILE x < 1.0D 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); + 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)/LONG(ln10); + END longlog; + + (* trigonometric functions; arguments in radians *) + + PROCEDURE sin(x: REAL): REAL; + BEGIN + RETURN SHORT(longsin(LONG(x))); + END sin; + + PROCEDURE longsin(x: LONGREAL): LONGREAL; + CONST + a2 = -0.1666666664D; + a4 = 0.0083333315D; + a6 = -0.0001984090D; + a8 = 0.0000027526D; + a10 = -0.0000000239D; + VAR + xsqr: LONGREAL; + neg: BOOLEAN; + BEGIN + neg := FALSE; + IF x < 0.0D THEN + neg := TRUE; + 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); + END; + IF x > LONG(halfpi) THEN + x := LONG(pi) - x; + 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; + END longsin; + + PROCEDURE cos(x: REAL): REAL; + BEGIN + RETURN SHORT(longcos(LONG(x))); + 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; + END longcos; + + PROCEDURE tan(x: REAL): REAL; + BEGIN + RETURN SHORT(longtan(LONG(x))); + END tan; + + PROCEDURE longtan(x: LONGREAL): LONGREAL; + VAR cosinus: LONGREAL; + BEGIN + cosinus := longcos(x); + IF cosinus = 0.0D THEN + (* ??? *) + RETURN 0.0D; + END; + RETURN longsin(x)/cosinus; + END longtan; + + PROCEDURE arcsin(x: REAL): REAL; + BEGIN + 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; + BEGIN + IF x < 0.0D THEN x := -x; END; + IF x > 1.0D THEN + (* ??? *) + RETURN 0.0D; + END; + RETURN LONG(halfpi) - + longsqrt(1.0D - x)*(a0+x*(a1+x*(a2+x*(a3+x*(a4+x*(a5+x*(a6+x*a7))))))); + END longarcsin; + + PROCEDURE arccos(x: REAL): REAL; + BEGIN + RETURN SHORT(longarccos(LONG(x))); + END arccos; + + PROCEDURE longarccos(x: LONGREAL): LONGREAL; + BEGIN + RETURN LONG(halfpi) - longarcsin(x); + END longarccos; + + PROCEDURE arctan(x: REAL): REAL; + BEGIN + RETURN SHORT(longarctan(LONG(x))); + END arctan; + + PROCEDURE longarctan(x: LONGREAL): LONGREAL; + 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; + VAR + xsqr: LONGREAL; + rev: BOOLEAN; + neg: BOOLEAN; + BEGIN + rev := FALSE; + neg := FALSE; + IF x < 0.0D THEN + neg := TRUE; + x := -x; + END; + IF x > 1.0D THEN + rev := TRUE; + x := 1.0D / x; + 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; + END; + IF neg THEN RETURN -x; END; + RETURN x; + END longarctan; + + (* hyperbolic functions *) + + 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 + (* ??? *) + RETURN 0.0D; + 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 + (* ??? *) + RETURN 0.0D; + 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; + +END Mathlib. diff --git a/lang/m2/libm2/RealConver.mod b/lang/m2/libm2/RealConver.mod index 8d510e175..9af9b9205 100644 --- a/lang/m2/libm2/RealConver.mod +++ b/lang/m2/libm2/RealConver.mod @@ -2,21 +2,30 @@ IMPLEMENTATION MODULE RealConversions; FROM FIFFEF IMPORT FIF, FEF; - PROCEDURE RealToString(r: REAL; + PROCEDURE RealToString(arg: REAL; + width, digits: INTEGER; + VAR str: ARRAY OF CHAR; + VAR ok: BOOLEAN); + BEGIN + LongRealToString(LONG(arg), width, digits, str, ok); + END RealToString; + + PROCEDURE LongRealToString(arg: LONGREAL; width, digits: INTEGER; VAR str: ARRAY OF CHAR; VAR ok: BOOLEAN); VAR pointpos: INTEGER; i: CARDINAL; ecvtflag: BOOLEAN; - intpart, fractpart: REAL; + r, intpart, fractpart: LONGREAL; ind1, ind2 : CARDINAL; sign: BOOLEAN; tmp : CHAR; ndigits: CARDINAL; - dummy, dig: REAL; + dummy, dig: LONGREAL; BEGIN + r := arg; DEC(width); IF digits < 0 THEN ecvtflag := TRUE; @@ -27,9 +36,10 @@ IMPLEMENTATION MODULE RealConversions; END; IF HIGH(str) < ndigits + 3 THEN str[0] := 0C; ok := FALSE; RETURN END; pointpos := 0; - sign := r < 0.0; + sign := r < 0.0D; IF sign THEN r := -r END; - r := FIF(r, 1.0, intpart); + r := FIF(r, 1.0D, intpart); + fractpart := r; pointpos := 0; ind1 := 0; ok := TRUE; @@ -37,9 +47,9 @@ IMPLEMENTATION MODULE RealConversions; Do integer part, which is now in "intpart". "r" now contains the fraction part. *) - IF intpart # 0.0 THEN + IF intpart # 0.0D THEN ind2 := 0; - WHILE intpart # 0.0 DO + WHILE intpart # 0.0D DO IF ind2 > HIGH(str) THEN IF NOT ecvtflag THEN str[0] := 0C; @@ -51,11 +61,11 @@ IMPLEMENTATION MODULE RealConversions; END; DEC(ind2); END; - dummy := FIF(FIF(intpart, 0.1, intpart),10.0, dig); - IF (dummy > 0.5) AND (dig < 9.0) THEN - dig := dig + 1.0; + dummy := FIF(FIF(intpart, 0.1D, intpart),10.0D, dig); + IF (dummy > 0.5D) AND (dig < 9.0D) THEN + dig := dig + 1.0D; END; - str[ind2] := CHR(TRUNC(dig+0.5) + ORD('0')); + str[ind2] := CHR(TRUNC(dig+0.5D) + ORD('0')); INC(ind2); INC(pointpos); END; @@ -70,10 +80,10 @@ IMPLEMENTATION MODULE RealConversions; END; ELSE INC(pointpos); - IF r > 0.0 THEN - WHILE r < 1.0 DO + IF r > 0.0D THEN + WHILE r < 1.0D DO fractpart := r; - r := r * 10.0; + r := r * 10.0D; DEC(pointpos); END; END; @@ -94,7 +104,7 @@ IMPLEMENTATION MODULE RealConversions; RETURN; END; WHILE ind1 <= ind2 DO - fractpart := FIF(fractpart, 10.0, r); + fractpart := FIF(fractpart, 10.0D, r); str[ind1] := CHR(TRUNC(r)+ORD('0')); INC(ind1); END; @@ -191,17 +201,26 @@ IMPLEMENTATION MODULE RealConversions; END; IF (ind1+1) <= HIGH(str) THEN str[ind1+1] := 0C; END; - END RealToString; + END LongRealToString; PROCEDURE StringToReal(str: ARRAY OF CHAR; VAR r: REAL; VAR ok: BOOLEAN); + VAR x: LONGREAL; + BEGIN + StringToLongReal(str, x, ok); + IF ok THEN + r := x; + END; + END StringToReal; - CONST BIG = 1.0E17; + PROCEDURE StringToLongReal(str: ARRAY OF CHAR; + VAR r: LONGREAL; VAR ok: BOOLEAN); + CONST BIG = 1.0D17; TYPE SETOFCHAR = SET OF CHAR; VAR pow10 : INTEGER; i : INTEGER; - e : REAL; + e : LONGREAL; ch : CHAR; signed: BOOLEAN; signedexp: BOOLEAN; @@ -209,11 +228,11 @@ IMPLEMENTATION MODULE RealConversions; PROCEDURE dig(ch: CARDINAL); BEGIN - IF r>BIG THEN INC(pow10) ELSE r:= 10.0*r + FLOAT(ch) END; + IF r>BIG THEN INC(pow10) ELSE r:= 10.0D*r + FLOATD(ch) END; END dig; BEGIN - r := 0.0; + r := 0.0D; pow10 := 0; iB := 0; ok := TRUE; @@ -276,10 +295,10 @@ IMPLEMENTATION MODULE RealConversions; pow10 := pow10 + i; END; IF pow10 < 0 THEN i := -pow10; ELSE i := pow10; END; - e := 1.0; + e := 1.0D; DEC(i); WHILE i >= 0 DO - e := e * 10.0; + e := e * 10.0D; DEC(i) END; IF pow10<0 THEN @@ -289,6 +308,6 @@ IMPLEMENTATION MODULE RealConversions; END; IF signed THEN r := -r; END; IF (iB <= HIGH(str)) AND (ORD(ch) > ORD(' ')) THEN ok := FALSE; END - END StringToReal; + END StringToLongReal; END RealConversions. -- 2.34.1