|
|
@@ -1,17 +1,27 @@
|
|
|
MODULE Reals;
|
|
|
|
|
|
-TYPE LONGREAL = REAL; (* REAL64 *) REAL = SHORTREAL; (* REAL32 *)
|
|
|
+IMPORT SYSTEM;
|
|
|
+
|
|
|
+TYPE
|
|
|
+ SHORTINT = SYSTEM.INT16;
|
|
|
+ INTEGER = SYSTEM.INT32;
|
|
|
+ HUGEINT = SYSTEM.INT64;
|
|
|
+ REAL = SYSTEM.REAL32;
|
|
|
+ CHAR = SYSTEM.CHAR16;
|
|
|
+
|
|
|
+ LONGREAL* = SYSTEM.REAL64;
|
|
|
|
|
|
VAR Done-: BOOLEAN;
|
|
|
|
|
|
-PROCEDURE Ten(e: INTEGER): REAL;
|
|
|
+PROCEDURE Ten(e: INTEGER): LONGREAL;
|
|
|
VAR r, power: LONGREAL;
|
|
|
-BEGIN r := 1.0; power := 10.0;
|
|
|
+BEGIN r := 1.0E0; power := 1.0E1;
|
|
|
WHILE e > 0 DO
|
|
|
- IF ODD(e) THEN r := r * power END;
|
|
|
- power := power * power; e := e DIV 2
|
|
|
- END ;
|
|
|
-RETURN SHORT(r) END Ten;
|
|
|
+ IF ODD(e) THEN r := r*power END;
|
|
|
+ power := power*power; e := SHORT(e DIV 2)
|
|
|
+ END;
|
|
|
+ RETURN r
|
|
|
+END Ten;
|
|
|
|
|
|
PROCEDURE LongVal*(IN s: ARRAY OF CHAR): LONGREAL;
|
|
|
VAR p, e: INTEGER; y, g: LONGREAL; neg, negE: BOOLEAN;
|
|
|
@@ -58,5 +68,297 @@ BEGIN
|
|
|
RETURN SHORT(LongVal(s))
|
|
|
END Val;
|
|
|
|
|
|
+(****)
|
|
|
+
|
|
|
+
|
|
|
+PROCEDURE digit(n: HUGEINT; VAR s: ARRAY OF CHAR; VAR i: INTEGER);
|
|
|
+BEGIN
|
|
|
+ DEC(i); s[i] := SHORT(CHR(n MOD 10 + 48));
|
|
|
+END digit;
|
|
|
+
|
|
|
+PROCEDURE Length(IN s: ARRAY OF CHAR): INTEGER;
|
|
|
+VAR n: INTEGER;
|
|
|
+BEGIN n := 0; WHILE (n < LEN(s)) & (s[n] # 0X) DO INC(n) END; RETURN n
|
|
|
+END Length;
|
|
|
+
|
|
|
+PROCEDURE prepend(IN t: ARRAY OF CHAR; VAR s: ARRAY OF CHAR; VAR i: INTEGER);
|
|
|
+ VAR j: INTEGER; l: INTEGER;
|
|
|
+BEGIN
|
|
|
+ l := Length(t); IF l > i THEN l := i END;
|
|
|
+ DEC(i, SHORT(l)); j := 0;
|
|
|
+ WHILE j < l DO s[i+j] := t[j]; INC(j) END
|
|
|
+END prepend;
|
|
|
+
|
|
|
+PROCEDURE -Entier64(x: LONGREAL): HUGEINT '(LONGINT)(x)';
|
|
|
+
|
|
|
+PROCEDURE RealP(x: LONGREAL; n: INTEGER; long: BOOLEAN; VAR s: ARRAY OF CHAR);
|
|
|
+(** Writes the long real number x to the end of the output
|
|
|
+ stream using an exponential form. If the textual representation of x
|
|
|
+ requires m characters (including a three-digit signed exponent), x is right
|
|
|
+ adjusted in a field of Max(n, m) characters padded with blanks at the left
|
|
|
+ end. A plus sign of the mantissa is not written. LONGREAL is 1/sign,
|
|
|
+ 11/exponent, 52/significand *)
|
|
|
+VAR
|
|
|
+ e: SHORTINT; (* Exponent field *)
|
|
|
+ f: HUGEINT; (* Fraction field *)
|
|
|
+ z: ARRAY 30 OF CHAR; (* Buffer built backwards *)
|
|
|
+ i: INTEGER; (* Index into z *)
|
|
|
+ el: INTEGER; (* Exponent length *)
|
|
|
+ x0: LONGREAL;
|
|
|
+ nn: BOOLEAN; (* Number negative *)
|
|
|
+ en: BOOLEAN; (* Exponent negative *)
|
|
|
+ m: HUGEINT; (* Mantissa digits *)
|
|
|
+ d: INTEGER; (* Significant digit count to display *)
|
|
|
+ dr: INTEGER; (* Number of insignificant digits that can be dropped *)
|
|
|
+ si: INTEGER; (* Runner pointing to s *)
|
|
|
+
|
|
|
+ PROCEDURE OutC(ch: CHAR);
|
|
|
+ BEGIN IF si < LEN(s) - 1 THEN s[si] := ch; INC(si) ELSE Done := FALSE END
|
|
|
+ END OutC;
|
|
|
+
|
|
|
+BEGIN si := 0; Done := TRUE;
|
|
|
+ e := SYSTEM.VAL(SHORTINT, (SYSTEM.VAL(HUGEINT, x) DIV 10000000000000L) MOD 800H);
|
|
|
+ f := SYSTEM.VAL(HUGEINT, x) MOD 10000000000000L;
|
|
|
+ nn := (SYSTEM.VAL(HUGEINT, x) < 0) & ~((e = 7FFH) & (f # 0)); (* Ignore sign on NaN *)
|
|
|
+ IF nn THEN DEC(n) END;
|
|
|
+
|
|
|
+ i := LEN(z);
|
|
|
+ IF e = 7FFH THEN (* NaN / Infinity *)
|
|
|
+ IF f = 0 THEN prepend('Infinity', z, i) ELSE prepend('NaN', z, i) END
|
|
|
+ ELSE
|
|
|
+ (* Calculate number of significant digits caller has proposed space for, and
|
|
|
+ number of digits to generate. *)
|
|
|
+ IF long THEN
|
|
|
+ el := 3;
|
|
|
+ dr := SHORT(n-6); (* Leave room for dp and '+D000' *)
|
|
|
+ IF dr > 17 THEN dr := 17 END; (* Limit to max useful significant digits *)
|
|
|
+ d := dr; (* Number of digits to generate *)
|
|
|
+ IF d < 15 THEN d := 15 END (* Generate enough digits to do trailing zero suppression *)
|
|
|
+ ELSE
|
|
|
+ el := 2;
|
|
|
+ dr := SHORT(n-5); (* Leave room for dp and '+E00' *)
|
|
|
+ IF dr > 9 THEN dr := 9 END; (* Limit to max useful significant digits *)
|
|
|
+ d := dr; (* Number of digits to generate *)
|
|
|
+ IF d < 6 THEN d := 6 END (* Generate enough digits to do trailing zero suppression *)
|
|
|
+ END;
|
|
|
+
|
|
|
+ IF e = 0 THEN
|
|
|
+ WHILE el > 0 DO DEC(i); z[i] := '0'; DEC(el) END;
|
|
|
+ DEC(i); z[i] := '+';
|
|
|
+ m := 0;
|
|
|
+ ELSE
|
|
|
+ IF nn THEN x := -x END;
|
|
|
+
|
|
|
+ (* Scale e to be an exponent of 10 rather than 2 *)
|
|
|
+ e := SHORT(SHORT(LONG(e - 1023) * 77 DIV 256));
|
|
|
+ IF e >= 0 THEN x := x / Ten(e) ELSE x := Ten(SHORT(-e)) * x END;
|
|
|
+ IF x >= 10.0E0 THEN x := 0.1E0 * x; INC(e) END;
|
|
|
+
|
|
|
+ (* Scale x to enough significant digits to reliably test for trailing
|
|
|
+ zeroes or to the amount of space available, if greater. *)
|
|
|
+ x0 := Ten(SHORT(d - 1));
|
|
|
+ x := x0 * x;
|
|
|
+ x := x + 0.5E0; (* Do not combine with previous line as doing so
|
|
|
+ introduces a least significant bit difference
|
|
|
+ between 32 bit and 64 bit builds. *)
|
|
|
+ IF x >= 10.0E0 * x0 THEN x := 0.1E0 * x; INC(e) END;
|
|
|
+
|
|
|
+ (* Generate the exponent digits *)
|
|
|
+ IF e < 0 THEN en := TRUE; e := SHORT(-e) ELSE en := FALSE END;
|
|
|
+ WHILE el > 0 DO digit(e, z, i); e := SHORT(e DIV 10); DEC(el) END;
|
|
|
+ DEC(i); IF en THEN z[i] := '-' ELSE z[i] := '+' END;
|
|
|
+
|
|
|
+ m := Entier64(x)
|
|
|
+ END;
|
|
|
+
|
|
|
+ DEC(i); IF long THEN z[i] := 'D' ELSE z[i] := 'E' END;
|
|
|
+
|
|
|
+ (* Drop trailing zeroes where caller proposes to use less space *)
|
|
|
+ IF dr < 2 THEN dr := 2 END;
|
|
|
+ WHILE (d > dr) & (m MOD 10 = 0) DO m := m DIV 10; DEC(d) END;
|
|
|
+
|
|
|
+ (* Render significant digits *)
|
|
|
+ WHILE d > 1 DO digit(m, z, i); m := m DIV 10; DEC(d) END;
|
|
|
+ DEC(i); z[i] := '.';
|
|
|
+ digit(m, z, i);
|
|
|
+ END;
|
|
|
+
|
|
|
+ (* Generate leading padding *)
|
|
|
+ DEC(n, SHORT(LEN(z)-i)); WHILE n > 0 DO OutC(' '); DEC(n) END;
|
|
|
+
|
|
|
+ (* Render prepared number from right end of buffer z *)
|
|
|
+ IF nn THEN OutC('-') END;
|
|
|
+ WHILE i < LEN(z) DO OutC(z[i]); INC(i) END;
|
|
|
+ s[si] := 0X
|
|
|
+END RealP;
|
|
|
+
|
|
|
+PROCEDURE Str*(x: REAL; n: INTEGER; VAR s: ARRAY OF CHAR);
|
|
|
+BEGIN RealP(x, n, FALSE, s)
|
|
|
+END Str;
|
|
|
+
|
|
|
+PROCEDURE LongStr*(x: LONGREAL; n: INTEGER; VAR s: ARRAY OF CHAR);
|
|
|
+BEGIN RealP(x, n, TRUE, s)
|
|
|
+END LongStr;
|
|
|
+
|
|
|
+(* Convert LONGREAL: Write positive integer value of x into array d.
|
|
|
+ The value is stored backwards, i.e. least significant digit
|
|
|
+ first. n digits are written, with trailing zeros fill.
|
|
|
+ On entry x has been scaled to the number of digits required. *)
|
|
|
+PROCEDURE ConvertL(x: LONGREAL; n: INTEGER; VAR d: ARRAY OF CHAR);
|
|
|
+ VAR i, j, k: HUGEINT;
|
|
|
+BEGIN
|
|
|
+ IF x < 0 THEN x := -x END;
|
|
|
+ k := 0;
|
|
|
+
|
|
|
+ IF (SIZE(INTEGER) < 8) & (n > 9) THEN
|
|
|
+ (* There are more decimal digits than can be held in a single INTEGER *)
|
|
|
+ i := ENTIER(x / 1000000000.0E0); (* The 10th and higher digits *)
|
|
|
+ j := ENTIER(x - (i * 1000000000.0E0)); (* The low 9 digits *)
|
|
|
+ (* First generate the low 9 digits. *)
|
|
|
+ IF j < 0 THEN j := 0 END;
|
|
|
+ WHILE k < 9 DO
|
|
|
+ d[k] := SHORT(CHR(j MOD 10 + 48)); j := j DIV 10; INC(k)
|
|
|
+ END;
|
|
|
+ (* Fall through to generate the upper digits *)
|
|
|
+ ELSE
|
|
|
+ (* We can generate all the digits in one go. *)
|
|
|
+ i := ENTIER(x);
|
|
|
+ END;
|
|
|
+
|
|
|
+ WHILE k < n DO
|
|
|
+ d[k] := SHORT(CHR(i MOD 10 + 48)); i := i DIV 10; INC(k)
|
|
|
+ END
|
|
|
+END ConvertL;
|
|
|
+
|
|
|
+PROCEDURE Expo(x: REAL): INTEGER;
|
|
|
+VAR i: INTEGER;
|
|
|
+BEGIN SYSTEM.GET(SYSTEM.ADR(x), i);
|
|
|
+ RETURN SHORT(ASH(i, -23) MOD 100H)
|
|
|
+END Expo;
|
|
|
+
|
|
|
+PROCEDURE LongExpo(x: LONGREAL): INTEGER;
|
|
|
+VAR i: HUGEINT;
|
|
|
+BEGIN SYSTEM.GET(SYSTEM.ADR(x), i);
|
|
|
+ RETURN SHORT(ASH(i, -52) MOD 800H)
|
|
|
+END LongExpo;
|
|
|
+
|
|
|
+PROCEDURE StrFix*(x: REAL; n, k: INTEGER; VAR s: ARRAY OF CHAR);
|
|
|
+CONST maxD = 9;
|
|
|
+VAR e, si, i, minus, p, N: INTEGER; x0: REAL;
|
|
|
+ d: ARRAY maxD OF CHAR;
|
|
|
+
|
|
|
+ PROCEDURE OutC(ch: CHAR);
|
|
|
+ BEGIN IF si < LEN(s) - 1 THEN s[si] := ch; INC(si) ELSE Done := FALSE END
|
|
|
+ END OutC;
|
|
|
+
|
|
|
+ PROCEDURE Chars(ch: CHAR; n: INTEGER);
|
|
|
+ BEGIN WHILE n > 0 DO OutC(ch); DEC(n) END
|
|
|
+ END Chars;
|
|
|
+
|
|
|
+ PROCEDURE DigitsOut(n: INTEGER);
|
|
|
+ BEGIN WHILE n > 0 DO DEC(i); OutC(d[i]); DEC(n) END
|
|
|
+ END DigitsOut;
|
|
|
+
|
|
|
+BEGIN Done := TRUE; si := 0; e := Expo(x);
|
|
|
+ IF k < 0 THEN k := 0 END;
|
|
|
+ IF e = 0 THEN
|
|
|
+ IF k = 0 THEN Chars(' ', n - 1); OutC('0')
|
|
|
+ ELSE Chars(' ', n - k - 2); OutC('0'); OutC('.'); Chars('0', k)
|
|
|
+ END
|
|
|
+ ELSIF e = 255 THEN Chars(' ', n - 3); OutC('N'); OutC('a'); OutC('N')
|
|
|
+ ELSE e := (e - 127) * 77 DIV 256;
|
|
|
+ IF x < 0 THEN minus := 1; x := -x ELSE minus := 0 END;
|
|
|
+ IF e >= 0 THEN (* x >= 1.0, 77/256 = log 2 *) x := SHORT(x / Ten(e))
|
|
|
+ ELSE (*x < 1.0*) x := SHORT(Ten(-e) * x)
|
|
|
+ END;
|
|
|
+ IF x >= 10.0 THEN x := 0.1 * x; INC(e) END;
|
|
|
+ (* 1 <= x < 10 *)
|
|
|
+ IF k + e >= maxD - 1 THEN k := maxD - 1 - e
|
|
|
+ ELSIF k + e < 0 THEN (* k := -e; !FIXME Why??? *) x := 0.0
|
|
|
+ END;
|
|
|
+ x0 := SHORT(Ten(k + e)); x := x0 * x + 0.5;
|
|
|
+ IF x >= 10.0 * x0 THEN INC(e) END;
|
|
|
+ (* e = no. of digits before decimal point *)
|
|
|
+ INC(e); i := k + e; ConvertL(x, i, d);
|
|
|
+ IF e > 0 THEN
|
|
|
+ IF k > 0 THEN p := 0 ELSE p := 1 END;
|
|
|
+ Chars(' ', n - e - k - 1 - minus + p);
|
|
|
+ IF minus # 0 THEN OutC('-') END;
|
|
|
+ DigitsOut(e);
|
|
|
+ IF k > 0 THEN OutC('.'); DigitsOut(k) END
|
|
|
+ ELSE
|
|
|
+ IF k + e > 0 THEN p := 0 ELSE p := 1 END;
|
|
|
+ N := n - k - 2 - minus + p;
|
|
|
+ IF (k > 0) & (k + e < 0) THEN DEC(N) END;
|
|
|
+ Chars(' ', N);
|
|
|
+ IF minus # 0 THEN OutC('-') END;
|
|
|
+ OutC('0');
|
|
|
+ IF k + e > 0 THEN OutC('.'); Chars('0', -e); DigitsOut(k + e)
|
|
|
+ ELSIF k > 0 THEN OutC('.'); Chars('0', k)
|
|
|
+ END
|
|
|
+ END
|
|
|
+ END;
|
|
|
+ s[si] := 0X
|
|
|
+END StrFix;
|
|
|
+
|
|
|
+PROCEDURE LongStrFix*(x: LONGREAL; n, k: INTEGER; VAR s: ARRAY OF CHAR);
|
|
|
+CONST maxD = 20;
|
|
|
+VAR e, si, i, minus, p, N: INTEGER; x0: LONGREAL;
|
|
|
+ d: ARRAY maxD OF CHAR;
|
|
|
+
|
|
|
+ PROCEDURE OutC(ch: CHAR);
|
|
|
+ BEGIN IF si < LEN(s) - 1 THEN s[si] := ch; INC(si) ELSE Done := FALSE END
|
|
|
+ END OutC;
|
|
|
+
|
|
|
+ PROCEDURE Chars(ch: CHAR; n: INTEGER);
|
|
|
+ BEGIN WHILE n > 0 DO OutC(ch); DEC(n) END
|
|
|
+ END Chars;
|
|
|
+
|
|
|
+ PROCEDURE DigitsOut(n: INTEGER);
|
|
|
+ BEGIN WHILE n > 0 DO DEC(i); OutC(d[i]); DEC(n) END
|
|
|
+ END DigitsOut;
|
|
|
+
|
|
|
+BEGIN Done := TRUE; si := 0; e := LongExpo(x);
|
|
|
+ IF k < 0 THEN k := 0 END;
|
|
|
+ IF e = 0 THEN
|
|
|
+ IF k = 0 THEN Chars(' ', n - 1); OutC('0')
|
|
|
+ ELSE Chars(' ', n - k - 2); OutC('0'); OutC('.'); Chars('0', k)
|
|
|
+ END
|
|
|
+ ELSIF e = 2047 THEN Chars(' ', n - 3); OutC('N'); OutC('a'); OutC('N')
|
|
|
+ ELSE e := (e - 1023) * 77 DIV 256;
|
|
|
+ IF x < 0 THEN minus := 1; x := -x ELSE minus := 0 END;
|
|
|
+ IF e >= 0 THEN (* x >= 1.0, 77/256 = log 2 *) x := x / Ten(e)
|
|
|
+ ELSE (*x < 1.0*) x := Ten(-e) * x
|
|
|
+ END;
|
|
|
+ IF x >= 10.0 THEN x := 0.1 * x; INC(e) END;
|
|
|
+ (* 1 <= x < 10 *)
|
|
|
+ IF k + e >= maxD - 1 THEN k := maxD - 1 - e
|
|
|
+ ELSIF k + e < 0 THEN (* k := -e; !FIXME Why??? *) x := 0.0
|
|
|
+ END;
|
|
|
+ x0 := Ten(k + e); x := x0 * x + 0.5;
|
|
|
+ IF x >= 10.0 * x0 THEN INC(e) END;
|
|
|
+ (* e = no. of digits before decimal point *)
|
|
|
+ INC(e); i := k + e; ConvertL(x, i, d);
|
|
|
+ IF e > 0 THEN
|
|
|
+ IF k > 0 THEN p := 0 ELSE p := 1 END;
|
|
|
+ Chars(' ', n - e - k - 1 - minus + p);
|
|
|
+ IF minus # 0 THEN OutC('-') END;
|
|
|
+ DigitsOut(e);
|
|
|
+ IF k > 0 THEN OutC('.'); DigitsOut(k) END
|
|
|
+ ELSE
|
|
|
+ IF k + e > 0 THEN p := 0 ELSE p := 1 END;
|
|
|
+ N := n - k - 2 - minus + p;
|
|
|
+ IF (k > 0) & (k + e < 0) THEN DEC(N) END;
|
|
|
+ Chars(' ', N);
|
|
|
+ IF minus # 0 THEN OutC('-') END;
|
|
|
+ OutC('0');
|
|
|
+ IF k + e > 0 THEN OutC('.'); Chars('0', -e); DigitsOut(k + e)
|
|
|
+ ELSIF k > 0 THEN OutC('.'); Chars('0', k)
|
|
|
+ END
|
|
|
+ END
|
|
|
+ END;
|
|
|
+ s[si] := 0X
|
|
|
+END LongStrFix;
|
|
|
+
|
|
|
BEGIN Done := TRUE
|
|
|
END Reals.
|
Arthur Yefimov