FIFFEF.def
InOut.def
Makefile
+Mathlib.def
MathLib0.def
Processes.def
RealInOut.def
+(*$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".
*)
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 ?
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
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:
DEFINITION MODULE MathLib0;
+(*
+ Exists for compatibility.
+ A more elaborate math lib can be found in Mathlib.def
+*)
PROCEDURE sqrt(x : REAL) : REAL;
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<sqrto2 THEN
- x := x + x;
- DEC(exp);
- END;
-
- z := (x-1.0)/(x+1.0);
- zsq := z*z;
-
- temp := ((p3*zsq + p2)*zsq + p1)*zsq + p0;
- temp := temp/(((zsq + q2)*zsq + q1)*zsq + q0);
- temp := temp*z + FLOAT(exp)*log2;
- RETURN temp;
+ RETURN Mathlib.ln(arg);
END ln;
-(*
- exp returns the exponential function of its
- floating-point argument.
-
- The coefficients are #1069 from Hart and Cheney. (22.35D)
-*)
-
-PROCEDURE floor(d: REAL): REAL;
-BEGIN
- IF d < 0.0 THEN
- d := -d;
- IF FIF(d, 1.0, d) # 0.0 THEN
- d := d + 1.0;
- END;
- d := -d;
- ELSE
- IF FIF(d, 1.0, d) # 0.0 THEN
- (* just ignore result of FIF *)
- ;
- END;
- END;
- RETURN d;
-END floor;
-
-PROCEDURE ldexp(fr: REAL; exp: INTEGER): REAL;
-VAR
- neg,i: INTEGER;
-BEGIN
- neg := 1;
- IF fr < 0.0 THEN
- fr := -fr;
- neg := -1;
- END;
- fr := FEF(fr, i);
- exp := exp + i;
- IF exp > 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;
--- /dev/null
+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.
--- /dev/null
+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.
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;
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;
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;
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;
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;
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;
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;
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;
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
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.
public / ack.git / commitdiff