EXPORTS Real =
BEGIN
PowTen: TYPE = RECORD [f: LONG CARDINAL, e: INTEGER];
TenTable: TYPE = RECORD [tens: ARRAY [0..13] OF PowTen, t26, t39: PowTen];
posTable: TenTable =
[tens:
[[20000000000B, 0], [24000000000B, 3], [31000000000B, 6], [37200000000B, 9],
[23420000000B, 13], [30324000000B, 16], [36411000000B, 19],
[23045500000B, 23], [27657020000B, 26], [35632624000B, 29],
[22500574400B, 33], [27220733500B, 36], [35065122420B, 39],
[22141163452B, 43]], t26: [24533722672B, 86], t39: [27405037645B, 129]];
negTable: TenTable =
[tens:
[[20000000000B, 0], [31463146315B, -4], [24365605075B, -7],
[20304467230B, -10], [32155613531B, -14], [24761326107B, -17],
[20615736406B, -20], [32657712326B, -24], [25363073422B, -27],
[21134057501B, -30], [33371577317B, -34], [25772777414B, -37],
[21457146011B, -40], [34113411502B, -44]], t26: [30604403045B, -87],
t39: [25616276613B, -130]];
iTenTable:
ARRAY [0..9]
OF
LONG
INTEGER =
[1, 10, 100, 1000, 10000, 100000, 1000000, 10000000, 100000000, 1000000000];
MulExtended:
PROC [x, y: IeeeInternal.Ext]
RETURNS [z: IeeeInternal.Ext] =
TRUSTED {
hi, lo: Basics.LongNumber;
z.exp ← x.exp + y.exp + 1;
z.det.sign ← x.det.sign # y.det.sign;
z.det.type ← normal;
z.det.sticky ← x.det.sticky OR y.det.sticky;
[hi, lo] ← IeeeInternal.Mul32[x.frac, y.frac]; -- normalize 64
WHILE
NOT IeeeInternal.BitOn[hi.highbits, IeeeInternal.HiBit]
DO
hi.lc ← hi.lc + hi.lc;
IF IeeeInternal.BitOn[lo.highbits, IeeeInternal.HiBit]
THEN
hi.lowbits ← PrincOpsUtils.BITOR[hi.lowbits, 1];
lo.lc ← lo.lc + lo.lc;
z.exp ← z.exp - 1;
ENDLOOP;
z.frac ← hi; -- Round to 32 bits.
IF lo.lc > IeeeInternal.HalfLC
OR
(lo.lc = IeeeInternal.HalfLC AND IeeeInternal.BitOn[z.frac.lowbits, 1]) THEN {
Overflow
z.frac.lc ← z.frac.lc + 1;
IF z.frac.lc < hi.lc
THEN {
z.frac.lc ← IeeeInternal.RShift1in1[z.frac.lc]; z.exp ← z.exp + 1; };
};
IF lo.lc # 0 THEN z.det.sticky ← TRUE;
};
RealToPair:
PUBLIC
PROC [
r: REAL, precision: CARDINAL ← Real.DefaultSinglePrecision]
RETURNS [type: Real.NumberType, fr: LONG INTEGER, exp10: INTEGER] = TRUSTED {
x, zz: IeeeInternal.Ext;
t1, t2, sign, normalized: BOOLEAN;
IeeeInternal.thisTimeExceptions ← Real.NoExceptions;
precision ← MAX[1, MIN[precision, 9]];
x ← IeeeInternal.Unpack[r];
sign ← x.det.sign;
SELECT x.det.type
FROM
nan => RETURN[nan, 0, 0];
infinity => {
fr ← IF sign THEN LAST[LONG INTEGER] ELSE FIRST[LONG INTEGER];
exp10 ← 99;
RETURN[infinity, fr, exp10];
};
zero => RETURN[zero, 0, 0];
ENDCASE => NULL;
x.det.sign ← FALSE;
normalized ← IeeeInternal.Normalized[x.frac.highbits];
exp10 ← IF normalized THEN (x.exp*3)/10 ELSE -38;
exp10 ← exp10 + 1 - precision;
DO
zz ← Scale[x, -exp10];
[v: fr, invalid: t1, overflow: t2] ← IeeeInternal.FixExtended[zz, rn];
IF t1 OR t2 THEN ERROR;
SELECT
TRUE
FROM
fr > iTenTable[precision] => exp10 ← exp10 + 1;
fr = iTenTable[precision] => {exp10 ← exp10 + 1; fr ← fr/10; EXIT; };
fr < iTenTable[precision - 1] =>
IF normalized THEN exp10 ← exp10 - 1 ELSE EXIT;
ENDCASE => EXIT;
ENDLOOP;
IF sign THEN fr ← -fr;
RETURN[normal, fr, exp10];
};
Scale:
PROC [x: IeeeInternal.Ext, exp10:
INTEGER]
RETURNS [y: IeeeInternal.Ext] =
TRUSTED {
table: TenTable;
mul: PowTen;
big: BOOLEAN;
IF exp10 = 0 THEN RETURN[x];
big ← exp10 < 0;
table ← IF big THEN negTable ELSE posTable;
exp10 ← ABS[exp10];
SELECT exp10
FROM
IN [1..13] => mul ← table.tens[exp10];
IN (13..26) => {
x ← MulExtended[x, CVExt[table.tens[13]]];
mul ← table.tens[exp10 - 13];
};
= 26 => mul ← table.t26;
IN (26..39] => {
x ← MulExtended[x, CVExt[table.t26]]; mul ← table.tens[exp10 - 26]; };
= 39 => mul ← table.t39;
IN (39..52] => {
x ← MulExtended[x, CVExt[table.t39]]; mul ← table.tens[exp10 - 39]; };
ENDCASE => {
WHILE exp10 > 52
DO
x ← MulExtended[x, CVExt[table.t39]]; exp10 ← exp10 - 39; ENDLOOP;
RETURN[Scale[x, IF big THEN -exp10 ELSE exp10]];
};
y ← CVExt[mul];
y ← MulExtended[x, y];
};
CVExt:
PROC [t: PowTen]
RETURNS [y: IeeeInternal.Ext] = {
y.det.sticky ← y.det.sign ← FALSE;
y.det.type ← normal;
y.frac.lc ← t.f;
y.exp ← t.e;
};
PairToReal:
PUBLIC
PROC [fr:
LONG
INTEGER, exp10:
INTEGER]
RETURNS [
REAL] =
TRUSTED {
y: IeeeInternal.Ext;
or: RealOps.RoundingMode ← IeeeInternal.fpmode.round;
IeeeInternal.fpmode.round ← rn;
IeeeInternal.thisTimeExceptions ← Real.NoExceptions;
IF fr = 0 THEN RETURN[Real.PlusZero];
y.exp ← 31;
y.det.sign ← fr < 0;
y.det.sticky ← FALSE;
y.det.type ← normal;
y.frac.li ← IF y.det.sign THEN -fr ELSE fr;
IeeeInternal.PostNormalize[@y];
y ← Scale[y, exp10];
IeeeInternal.StepTwo[@y];
IeeeInternal.fpmode.round ← or;
RETURN[IeeeInternal.Pack[@y]];
};
END.