IeeeIOA.mesa - Utilities for IO
Copyright © 1985 by Xerox Corporation. All rights reserved.
Rovner On May 4, 1983 9:59 am
Stewart, January 14, 1984 4:18 pm
Levin, August 8, 1983 4:38 pm
Russ Atkinson (RRA) May 23, 1985 0:56:38 am PDT
DIRECTORY
Basics USING [LongNumber],
Ieee USING [BitOn, Ext, HalfLC, HiBit, Mul32, Pack, PostNormalize, RShift1in1, SingleReal, StepTwo],
PrincOpsUtils USING [BITOR],
Real USING [MaxSinglePrecision, NumberType, RoundLI],
RealExceptions USING [ClearExceptions],
RealOps USING [DefMode, Mode, SetMode];
IeeeIOA:
CEDAR
PROGRAM
IMPORTS Ieee, PrincOpsUtils, Real, RealExceptions, RealOps
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]];
MulExtended:
PROC [x, y: Ieee.Ext]
RETURNS [z: Ieee.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] ← Ieee.Mul32[x.frac, y.frac]; -- normalize 64
WHILE
NOT Ieee.BitOn[hi.highbits, Ieee.HiBit]
DO
hi.lc ← hi.lc + hi.lc;
IF Ieee.BitOn[lo.highbits, Ieee.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 > Ieee.HalfLC
OR
(lo.lc = Ieee.HalfLC AND Ieee.BitOn[z.frac.lowbits, 1]) THEN {
Overflow
z.frac.lc ← z.frac.lc + 1;
IF z.frac.lc < hi.lc
THEN {
z.frac ← Ieee.RShift1in1[z.frac];
z.exp ← z.exp + 1; };
};
IF lo.lc # 0 THEN z.det.sticky ← TRUE;
};
myRealToPairTable:
REF RealTable ← NewRealTable[];
RealTableIndex: TYPE = [-38..38];
RealTable: TYPE = ARRAY RealTableIndex OF REAL;
myTensTable: TensTableRef ←
NEW[TensTableArray ← [
1, 10, 100, 1000, 10000, 100000, 1000000, 10000000, 100000000, 1000000000]];
TensTableRef: TYPE = REF TensTableArray;
TensTableArray: TYPE = ARRAY [0..Real.MaxSinglePrecision] OF INT;
NewRealTable:
PROC
RETURNS [ref:
REF RealTable] = {
ref ← NEW[RealTable];
FOR i: RealTableIndex
IN RealTableIndex
DO
ref[i] ← PairToReal[1, i];
ENDLOOP;
};
RealToPair:
PUBLIC
PROC [r:
REAL, precision:
NAT]
RETURNS [type: Real.NumberType, fr:
INT, exp10:
INTEGER] = {
lo: RealTableIndex ← -38;
hi: RealTableIndex ← LAST[RealTableIndex];
abs: REAL ← ABS[r];
single: Ieee.SingleReal ← LOOPHOLE[abs, Ieee.SingleReal];
adjust: INTEGER ← 0;
sexp: [0..377B] ← single.exp;
exp2: INTEGER ← sexp-200B;
type ← normal;
SELECT sexp
FROM
< 40B => {
SELECT
TRUE
FROM
r = 0.0 =>
0.0 is a special case.
RETURN [zero, 0, 0];
ENDCASE => {
Denormalized or VERY small numbers
floor: REAL ← myRealToPairTable[precision-38];
WHILE abs < floor
DO
abs ← abs * 1e10;
adjust ← adjust - 10;
ENDLOOP;
hi ← -18;
};
};
>= 340B => {
IF sexp # 377B
THEN
Normal (but but quite large) numbers. We need to watch out for edge effects in estimating lo and hi bounds for the search.
lo ← exp2*3/10
ELSE {
Infinity or Not a Number
IF r < 0 THEN fr ← FIRST[INT] ELSE fr ← LAST[INT];
exp10 ← 99;
IF single.m1 = 0
AND single.m2 = 0
THEN type ← infinity
ELSE type ← nan;
RETURN;
};
};
ENDCASE => {
Normal numbers, easy to estimate the lo and hi bounds.
lo ← exp2*3/10-2;
hi ← lo + 4;
};
Take precision into the range [1..Real.MaxSinglePrecision]. We always give back at least one digit of precision.
SELECT precision
FROM
0 => precision ← 1;
> Real.MaxSinglePrecision => precision ← Real.MaxSinglePrecision;
ENDCASE;
Perform a binary search for the best power of ten.
DO
mid: RealTableIndex ← (lo+hi)/2;
guess: REAL ← myRealToPairTable[mid];
IF mid = lo
OR mid = hi
THEN {
The binary chop has converged. The correct value of mid is either lo or hi. So, a simple test should suffice.
IF mid = hi THEN mid ← lo;
IF mid < hi
THEN {
Well, we might be too low by one index.
next: REAL ← myRealToPairTable[mid+1];
IF next < abs THEN mid ← mid + 1;
};
mid ← mid-precision+1;
We believe that mid IN RealTableIndex at this point, because of the check for adjusting denormalized numbers up within the floor. If this is not correct, we deserve to get a bounds fault.
exp10 ← mid+adjust;
abs ← abs / myRealToPairTable[mid];
EXIT;
};
IF guess <= abs THEN {lo ← mid; LOOP};
hi ← mid;
ENDLOOP;
fr ← Real.RoundLI[abs];
IF fr = myTensTable[precision]
THEN {
Rats! The Real.RoundLI operation rounded up to bump the # of precision, so we have to adjust the # of precision down, and adjust the exponent as well.
fr ← myTensTable[precision-1];
exp10 ← exp10 + 1;
};
IF r < 0.0 THEN fr ← - fr;
};
Scale:
PROC [x: Ieee.Ext, exp10:
INTEGER]
RETURNS [y: Ieee.Ext] =
TRUSTED {
table: TenTable;
mul: PowTen;
big: BOOL;
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: Ieee.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:
INT, exp10:
INTEGER]
RETURNS [r:
REAL] =
TRUSTED {
RealExceptions.ClearExceptions[];
IF fr = 0
THEN r ← 0.0
ELSE {
y: Ieee.Ext;
old: RealOps.Mode ← RealOps.SetMode[RealOps.DefMode];
{
ENABLE UNWIND => [] ← RealOps.SetMode[old];
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;
Ieee.PostNormalize[@y];
y ← Scale[y, exp10];
Ieee.StepTwo[@y];
r ← Ieee.Pack[@y];
};
[] ← RealOps.SetMode[old];
};
};
August25, 1980 4:28 PM, LStewart; formatting and qualification
August 27, 1982 1:04 pm, L. Stewart, CEDAR
January 14, 1984 4:53 pm, L. Stewart, changeTo Ieee
Russ Atkinson (RRA) February 19, 1985 5:05:51 pm PST
fixed RealToPair to be more accurate and faster