-- IeeePack.mesa
-- Ieee Floating Point Package for compiler (Ieee + IeeeIOA + IeeeUtil)
-- Last Modified by Satterthwaite, May 31, 1982 12:04 pm
DIRECTORY
Inline: TYPE USING [BITAND, BITOR, BITSHIFT, BITXOR, LongMult, LongNumber],
Real: FROM "IeeeFloat" USING [
Exception, ExceptionFlags, Extended, MinusInfinity, MinusZero, Mode, NoExceptions,
NumberType, PlusInfinity, PlusZero, TrapNonTrappingNaN, TrapTrappingNaN];
IeeePack: PROGRAM
IMPORTS Inline
EXPORTS Real = {
OPEN Inline, Real;
RealException: PUBLIC SIGNAL [flags: ExceptionFlags] = CODE;
-- from Ieee.Mesa
Details: TYPE = MACHINE DEPENDENT RECORD [
sign: BOOL,
sticky: BOOL,
blank: [0..7777b],
type: Real.NumberType];
Ext: TYPE = RECORD [
det: Details,
exp: INTEGER,
frac: Inline.LongNumber];
-- constants
HiBit: CARDINAL = 100000b;
ExponentBias: INTEGER = 127;
ExponentMask: INTEGER = 077600b;
ExponentShift: INTEGER = 7;
HiFractionMask: INTEGER = 177b;
FractionShift: INTEGER = 8;
LeastSignificandBit: INTEGER = 0400b;
HiddenBit: CARDINAL = 100000b;
StickyBits: CARDINAL = 77b;
HalfLC: LONG CARDINAL = 20000000000b;
SignBit: CARDINAL = 100000b;
ExpSingleMax: INTEGER = 127;
ExpSingleMin: INTEGER = -126;
DenormalizedExponent: INTEGER = -127;
NaNExponent: INTEGER = 128;
BiasAdjust: INTEGER = 192;
-- global variables
fpmode: Mode = [nm: warning, round: rn];
enableFlags: ExceptionFlags = [invalidOperation: TRUE, divisionByZero: TRUE, overflow: TRUE];
thisTimeExceptions: ExceptionFlags;
LN: PROC [r: LONG UNSPECIFIED] RETURNS [LongNumber] = INLINE {
RETURN[LOOPHOLE[r, LongNumber]]};
BitOn: PROC [a, b: UNSPECIFIED] RETURNS [BOOL] = INLINE {
RETURN[BITAND[a, b]#0]};
Normalized: PROC [g: INTEGER] RETURNS [BOOL] = INLINE {
RETURN [BITAND[g, HiddenBit]#0]};
-- Add with Carry
ADC3: PROC [a, b, c: CARDINAL] RETURNS [CARDINAL, CARDINAL] = INLINE {
s: LongNumber;
s.lc ← LONG[a]+LONG[b]+LONG[c];
RETURN[s.highbits, s.lowbits]};
-- from IeeeIOA
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: Ext] RETURNS [z: Ext] = {
hi, lo: 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] ← Mul32[x.frac, y.frac];
-- normalize 64
WHILE NOT BitOn[hi.highbits, HiBit] DO
hi.lc ← hi.lc+hi.lc;
IF BitOn[lo.highbits, HiBit] THEN hi.lowbits ← 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>HalfLC OR (lo.lc=HalfLC AND BitOn[z.frac.lowbits, 1]) THEN {
-- Overflow
z.frac.lc ← z.frac.lc+1;
IF z.frac.lc<hi.lc THEN {
z.frac.lc ← RShift1in1[z.frac.lc];
z.exp ← z.exp+1};
};
IF lo.lc#0 THEN z.det.sticky ← TRUE};
Scale: PROC [x: Ext, exp10: INTEGER] RETURNS [y: Ext] = {
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: 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] = {
y: Ext;
-- or: RoundingMode ← fpmode.round;
-- fpmode.round ← rn;
thisTimeExceptions ← NoExceptions;
IF fr=0 THEN RETURN[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;
PostNormalize[@y];
y ← Scale[y, exp10];
StepTwo[@y];
-- fpmode.round ← or;
RETURN[Pack[@y]]};
-- from IeeeFloatA (adapted)
RealToExtended: PUBLIC PROC [a: REAL] RETURNS [Real.Extended] = {
ext: Ext = Unpack[a];
RETURN [[type: ext.det.type, sign: ext.det.sign, exp: ext.exp, frac: ext.frac.lc]]};
ExtendedToReal: PUBLIC PROC [e: Real.Extended] RETURNS [REAL] = {
ext: Ext ← [
det: [sign: e.sign, sticky: FALSE, blank: 0, type: e.type],
exp: e.exp, frac: [lc[e.frac]]];
RETURN [Pack[@ext]]};
Negate: PUBLIC PROC [a: REAL] RETURNS [REAL] = {
IF LN[a] = LN[PlusZero] THEN RETURN [PlusZero];
-- Negate a
LOOPHOLE[a, LongNumber].highbits ← BITXOR[LN[a].highbits, SignBit];
RETURN [a]};
Abs: PUBLIC PROC [a: REAL] RETURNS [REAL] = {
RETURN [IF BitOn[LN[a].highbits, SignBit] THEN Negate[a] ELSE a]};
-- from IEEEUtil
SetInexactResult: PROC = {
thisTimeExceptions[inexactResult] ← TRUE};
SetInvalidOperation: PROC = {
thisTimeExceptions[invalidOperation] ← TRUE};
SetUnderflow: PROC [z: POINTER TO Ext] = {
thisTimeExceptions[underflow] ← TRUE};
SetOverflow: PROC [z: POINTER TO Ext] = {
thisTimeExceptions[overflow] ← TRUE};
-- Separate the packed REAL into its component elements
Unpack: PUBLIC PROC [r: REAL] RETURNS [z: Ext] = {
z.det.sticky ← FALSE;
z.det.sign ← BitOn[LN[r].highbits, HiBit]; -- first r.j. the exponent
z.exp ← Inline.BITSHIFT[
Inline.BITAND[LN[r].highbits, ExponentMask], -ExponentShift];
z.exp ← z.exp - ExponentBias;
z.det.type ← normal;
z.frac.li ← LN[r].li;
z.frac.highbits ← Inline.BITAND[HiFractionMask, z.frac.highbits];
SELECT z.exp FROM
= DenormalizedExponent => {
-- denormalized or zero
IF z.frac.li = 0 THEN z.det.type ← zero
ELSE {
z.exp ← ExpSingleMin;
z.frac.li ← LongShift[z.frac.li, FractionShift];
IF fpmode.nm = normalizing THEN PostNormalize[@z]};
};
= NaNExponent => {
-- infinity or nan
IF z.frac.li = 0 THEN z.det.type ← infinity
ELSE {
z.det.type ← nan;
IF z.frac.lc = Real.TrapTrappingNaN THEN SetInvalidOperation[]};
};
ENDCASE => {
z.frac.li ← LongShift[z.frac.li, FractionShift];
z.frac.highbits ← Inline.BITOR[HiddenBit, z.frac.highbits]};
};
-- Stuff the components back into the packed REAL.
Pack: PROC [z: POINTER TO Ext] RETURNS [r: REAL] = {
trap: BOOL ← FALSE;
i: Exception;
IF thisTimeExceptions#NoExceptions THEN {
-- Possible typo in standard here!!!
IF thisTimeExceptions[invalidOperation] THEN {
z.det.type ← nan;
IF z.frac.lc=0 THEN z.frac.lc ← TrapNonTrappingNaN};
FOR i IN Exception DO
trap ← trap OR (enableFlags[i] AND thisTimeExceptions[i]);
ENDLOOP ;
IF trap THEN
SIGNAL RealException[flags: thisTimeExceptions];
FixupProcedure[z]};
RETURN[UsualPack[z↑]]};
UsualPack: PROC [z: Ext] RETURNS [REAL] = {
SELECT z.det.type FROM
zero => RETURN [IF z.det.sign THEN MinusZero ELSE PlusZero];
infinity => RETURN [IF z.det.sign THEN MinusInfinity ELSE PlusInfinity];
nan => z.exp ← NaNExponent;
ENDCASE => z.frac.li ← LongShift[z.frac.li, -FractionShift];
-- clear hidden bit
z.frac.highbits ← BITAND[z.frac.highbits, HiFractionMask];
IF z.exp NOT IN [-127..128] THEN ERROR;
z.exp ← BITSHIFT[z.exp+ExponentBias, ExponentShift];
z.frac.highbits ← BITOR[z.frac.highbits, z.exp];
IF z.det.sign THEN z.frac.highbits ← BITOR[z.frac.highbits, SignBit];
RETURN [LOOPHOLE[z.frac.li, REAL]]};
FixupProcedure: PROC [vp: POINTER TO Ext] = {
IF thisTimeExceptions[underflow] THEN {
DeNormalize[vp, vp.exp-ExpSingleMin];
vp.exp ← DenormalizedExponent};
IF thisTimeExceptions[overflow] THEN
vp.det.type ← infinity;
};
Round: PROC [z: POINTER TO Ext] = {
temp: LONG CARDINAL;
-- A 32 bit extended is considered to have its hidden bit as bit 0, significand as the next 29 bits, then G, then R, then S as a separate boolean. GRS computs the 3 bit number formed by the concatenation of G, R, and IF S THEN 1 ELSE 0.
GRS: PROC [g: INTEGER] RETURNS [INTEGER] = {
-- StickyBits are the set used for rounding to 24 bits
s: BOOL ← BitOn[g, StickyBits] OR z.det.sticky;
-- Should flush these constants
g ← BITSHIFT[BITAND[g, 300b], -5];
IF s THEN g ← g+1;
RETURN[g]};
-- The fraction should be normalized here!
grs: INTEGER ← GRS[z.frac.lowbits];
IF grs=0 THEN RETURN;
SetInexactResult[];
BEGIN
SELECT fpmode.round FROM
rn => {
IF z.det.sign THEN grs ← 8-grs;
IF grs>4 OR (grs=4 AND BitOn[z.frac.lowbits, LeastSignificandBit]) THEN
GOTO Plus1;
};
rz => NULL;
rm => IF z.det.sign THEN GOTO Plus1;
rp => IF NOT z.det.sign THEN GOTO Plus1;
ENDCASE => ERROR;
EXITS
Plus1 => {
temp ← z.frac.lc;
z.frac.lc ← z.frac.lc+LeastSignificandBit;
IF z.frac.lc<=temp THEN { -- better overflow detect!
z.frac.li ← RShift1in1[z.frac.li];
z.exp ← z.exp+1};
};
END};
RShift1in1: PROC [z: LONG UNSPECIFIED] RETURNS [LONG UNSPECIFIED] = {
vl: LongNumber;
vl.lc ← z;
vl.lowbits ← BITSHIFT[vl.lowbits, -1];
IF BitOn[vl.highbits, 1] THEN vl.lowbits ← BITOR[vl.lowbits, HiBit];
vl.highbits ← BITOR[HiddenBit, BITSHIFT[vl.highbits, -1]];
RETURN [vl.lc]};
LShift: PROC [z: LONG UNSPECIFIED] RETURNS [LONG UNSPECIFIED] = {
vl: LongNumber;
vl.lc ← z;
vl.highbits ← BITSHIFT[vl.highbits, 1];
IF BitOn[vl.lowbits, HiBit] THEN vl.highbits ← BITOR[vl.highbits, 1];
vl.lowbits ← BITSHIFT[vl.lowbits, 1];
RETURN [vl.lc]};
StepTwo: PROC [z: POINTER TO Ext] = {
IF z.det.type#normal THEN RETURN;
IF z.exp<=ExpSingleMin THEN {
IF z.exp<ExpSingleMin OR (NOT Normalized[z.frac.highbits]) THEN SetUnderflow[z]};
Round[z];
IF (NOT Normalized[z.frac.highbits]) AND z.exp#DenormalizedExponent AND (NOT thisTimeExceptions[underflow]) THEN SetInvalidOperation[]
ELSE IF z.exp>ExpSingleMax THEN SetOverflow[z]};
Mul32: PROC [x, y: LongNumber] RETURNS [LongNumber, LongNumber] = {
hi, lo, t1, t2: LongNumber;
cy: CARDINAL;
lo.lc ← LongMult[x.lowbits, y.lowbits];
hi.lc ← LongMult[x.highbits, y.highbits];
t1.lc ← LongMult[x.highbits, y.lowbits];
t2.lc ← LongMult[x.lowbits, y.highbits];
[cy, lo.highbits] ← ADC3[lo.highbits, t1.lowbits, t2.lowbits];
hi.lc ← hi.lc+t1.highbits+t2.highbits+cy;
RETURN [hi, lo]};
-- Post Normalize. S does not participate
PostNormalize: PROC [z: POINTER TO Ext] = {
IF z.frac.lc=0 THEN ERROR;
WHILE BITAND[z.frac.highbits, HiddenBit]=0 DO
z.frac.lc ← LShift[z.frac.lc];
z.exp ← z.exp-1;
ENDLOOP;
};
-- positive count is left shift, negative is right shift
LongShift: PROC [z: LONG UNSPECIFIED, count: INTEGER] RETURNS [LONG UNSPECIFIED] = {
xl, vl: LongNumber;
vl.lc ← z;
SELECT count FROM
NOT IN (-32..32) => RETURN [0];
0 => RETURN [vl.lc];
IN [16..32) => {
xl.highbits ← BITSHIFT[vl.lowbits,count-16];
xl.lowbits ← 0};
IN (0..16) => {
xl.highbits ← BITSHIFT[vl.highbits, count] + BITSHIFT[vl.lowbits,count-16];
xl.lowbits ← BITSHIFT[vl.lowbits, count]};
IN (-16..0) => {
xl.highbits ← BITSHIFT[vl.highbits, count];
xl.lowbits ← BITSHIFT[vl.lowbits, count] + BITSHIFT[vl.highbits,count+16]};
IN (-32..-16] => {
xl.highbits ← 0;
xl.lowbits ← BITSHIFT[vl.highbits,count+16]};
ENDCASE => ERROR;
RETURN[xl.lc]};
-- DeNormalize is much like LongShift, except that it maintains the sticky bits on the right. And it only shifts right.
DeNormalize: PROC [z: POINTER TO Ext, count: INTEGER] = {
sMask: ARRAY [1..16) OF CARDINAL = [
1b, 3b, 7b, 17b, 37b, 77b, 177b, 377b, 777b, 1777b, 3777b, 7777b, 17777b, 37777b, 77777b];
-- Mask off everything but the bits contributing to S, then set or clear S in the result.
IF count>0 THEN ERROR;
IF count=0 THEN RETURN;
count ← -count;
z.det.sticky ← SELECT count FROM
=16 => z.frac.lowbits#0,
IN [1..16) => BitOn[z.frac.lowbits, sMask[count]],
IN (16..32) => z.frac.lowbits#0 OR BitOn[z.frac.highbits, sMask[count-16]],
ENDCASE => z.frac.lc#0;
z.frac.lc ← LongShift[z.frac.lc, -count]};
}.