-- IeeeFloatA.mesa
-- Last Modified: Stewart 17-Nov-81 1:14:34
-- Mesa implementation of floating point ops +, -, *, /, Rem
DIRECTORY
Ieee,
Inline,
Real USING [
AddInfinityNaN, DivideInfinityNaN, MinusZero, MultiplyInfinityNaN,
NoExceptions, NumberType, PlusZero],
RealOps;
IeeeFloatA: PROGRAM IMPORTS Ieee, Inline EXPORTS RealOps =
BEGIN OPEN Ieee, Inline, Real, RealOps;
AddNormal: PROC [x, y: POINTER TO Ext] = {
ediff: INTEGER ← x.exp - y.exp;
normalized: BOOLEAN;
IF ediff > 0 THEN DeNormalize[y, -ediff]
ELSE {DeNormalize[x, ediff]; x.exp ← y.exp; };
normalized ← Normalized[BITOR[x.frac.highbits, y.frac.highbits]];
IF x.det.sign = y.det.sign THEN {
cy: CARDINAL;
[cy, x.frac.lowbits] ← ADC2[x.frac.lowbits, y.frac.lowbits];
[cy, x.frac.highbits] ← ADC3[x.frac.highbits, y.frac.highbits, cy];
IF cy # 0 THEN {
IF BitOn[x.frac.lowbits, 1] THEN
x.frac.lowbits ← Inline.BITOR[x.frac.lowbits, 2];
x.frac.li ← RShift1in1[x.frac.li];
x.exp ← x.exp + 1;
};
}
ELSE {
IF x.frac.lc >= y.frac.lc THEN x.frac.lc ← x.frac.lc - y.frac.lc
ELSE {x.frac.lc ← y.frac.lc - x.frac.lc; x.det.sign ← y.det.sign; };
IF x.frac.highbits = 0 AND NOT BitOn[x.frac.lowbits, LowSignificandMask]
THEN {x.det.sign ← fpmode.round = rm; x.det.type ← zero; }
ELSE {IF normalized THEN PostNormalize[x]; };
};
IF y.det.sticky THEN x.det.sticky ← TRUE;
};
FAdd: PUBLIC SAFE PROC [a, b: REAL, m: Mode] RETURNS [REAL] = TRUSTED {
x, y: Ext;
{
fpmode ← m;
thisTimeExceptions ← NoExceptions;
x ← Unpack[a];
y ← Unpack[b];
IF BitOn[funny, BITOR[x.det, y.det]] THEN
SELECT TRUE FROM
x.det.type = nan => GOTO RetX;
y.det.type = nan => GOTO RetY;
x.det.type # infinity AND y.det.type # infinity => NULL;
x.det.type = infinity AND y.det.type = infinity =>
IF fpmode.im = affine AND x.det.sign = y.det.sign THEN GOTO RetA
ELSE {
SetInvalidOperation[];
x.det.type ← nan;
x.frac.lc ← AddInfinityNaN;
GOTO RetX;
};
x.det.type = infinity => GOTO RetA;
y.det.type = infinity => GOTO RetB;
ENDCASE => ERROR
ELSE
SELECT NormalType[x.det, y.det] FROM
TNnn => NULL;
TNnz => GOTO RetA;
TNzn => GOTO RetB;
TNzz =>
RETURN[
IF fpmode.round = rm AND x.det.sign AND y.det.sign THEN MinusZero
ELSE PlusZero];
ENDCASE => ERROR;
AddNormal[@x, @y];
StepTwo[@x];
GOTO RetX;
EXITS
RetA => RETURN[a];
RetB => RETURN[b];
RetX => RETURN[Pack[@x]];
RetY => RETURN[Pack[@y]];
};
};
FSub: PUBLIC SAFE PROC [a, b: REAL, m: Mode] RETURNS [REAL] = TRUSTED {
-- Negate b
LOOPHOLE[b, LongNumber].highbits ← BITXOR[LN[b].highbits, SignBit];
RETURN[FAdd[a, b, m]];
};
FMul: PUBLIC SAFE PROC [a, b: REAL, m: Mode] RETURNS [REAL] = TRUSTED {
x, y: Ext;
lo: LongNumber;
{
fpmode ← m;
thisTimeExceptions ← NoExceptions;
x ← Unpack[a];
y ← Unpack[b];
x.det.sign ← x.det.sign # y.det.sign; -- due to different bias
x.exp ← x.exp + y.exp + 1;
IF BitOn[funny, BITOR[x.det, y.det]] THEN
SELECT TRUE FROM
x.det.type = nan => GOTO RetX;
y.det.type = nan => GOTO RetY;
x.det.type = infinity OR y.det.type = infinity => {
IF x.det.type = zero OR y.det.type = zero THEN {
SetInvalidOperation[];
x.det.type ← nan;
x.frac.lc ← MultiplyInfinityNaN;
}; -- Possible test for unnormal zero later.
x.det.type ← infinity;
GOTO RetX;
};
ENDCASE => NULL -- one or more denormalized
ELSE
SELECT NormalType[x.det, y.det] FROM
TNnn => NULL;
TNnz => GOTO RetZero;
TNzn => GOTO RetZero;
TNzz => GOTO RetZero;
ENDCASE => ERROR;
[x.frac, lo] ← Mul32[x.frac, y.frac]; -- normalize (double normalize??)
IF x.frac.lc=0 THEN {
x.frac.lc ← lo.lc;
x.exp ← x.exp-32;
}
ELSE {IF lo.lc # 0 THEN x.det.sticky ← TRUE;};
PostNormalize[@x];
StepTwo[@x];
GOTO RetX;
EXITS
RetX => RETURN[Pack[@x]];
RetY => RETURN[Pack[@y]];
RetZero => RETURN[IF x.det.sign THEN MinusZero ELSE PlusZero];
};
};
FDiv: PUBLIC SAFE PROC [a, b: REAL, m: Mode] RETURNS [c: REAL] = TRUSTED {
x, y: Ext;
unnormalized: BOOLEAN ← FALSE;
{
fpmode ← m;
thisTimeExceptions ← NoExceptions;
x ← Unpack[a];
y ← Unpack[b];
x.det.sign ← x.det.sign # y.det.sign;
{
IF BitOn[funny, BITOR[x.det, y.det]] THEN {
SELECT TRUE FROM
x.det.type = nan => GOTO RetX;
y.det.type = nan => GOTO RetY;
x.det.type = infinity AND y.det.type = infinity => GOTO ZeroDivide;
-- possible test for unnormal zero on this one
y.det.type = infinity => GOTO RetZero;
x.det.type = infinity => GOTO RetX;
ENDCASE => IF NOT Normalized[y.frac.highbits] THEN GOTO ZeroDivide;
}
ELSE
SELECT NormalType[x.det, y.det] FROM
TNnn => IF Normalized[y.frac.highbits] THEN NULL ELSE GOTO ZeroDivide;
-- this is usual flow
TNnz => GOTO RetInf;
TNzn => GOTO RetZero;
TNzz => GOTO ZeroDivide;
ENDCASE => ERROR;
EXITS
ZeroDivide => {
SetInvalidOperation[];
x.det.type ← nan;
x.frac.lc ← DivideInfinityNaN;
GOTO RetX;
};
RetInf => {SetDivisionByZero[]; x.det.type ← infinity; GOTO RetX; };
};
-- Have to make sure the numbers are normalized for this to work, so we remember if either was unnormalized (warning).
IF NOT Normalized[x.frac.highbits] THEN {
PostNormalize[@x]; unnormalized ← TRUE; };
x ← Divide[x, y];
StepTwo[@x];
GOTO RetX;
EXITS
RetZero => RETURN[IF x.det.sign THEN MinusZero ELSE PlusZero];
RetX => RETURN[Pack[@x]];
RetY => RETURN[Pack[@y]];
};
};
Divide: PROC [a, b: Ext] RETURNS [y: Ext] = {
cy: BOOLEAN;
Step: PROC = INLINE {
y.frac.lc ← y.frac.lc + y.frac.lc;
IF cy OR a.frac.lc >= b.frac.lc THEN {
a.frac.lc ← a.frac.lc - b.frac.lc;
y.frac.lowbits ← y.frac.lowbits + 1; -- can't carry!
};
cy ← BitOn[a.frac.highbits, HiBit];
a.frac.lc ← a.frac.lc + a.frac.lc;
}; -- initialize
IF b.frac.lc = 0 THEN ERROR;
y ← a;
y.det ←
[sign: a.det.sign, sticky: FALSE, blank: 0, type: Real.NumberType[normal]];
y.exp ← a.exp - b.exp;
y.frac.lc ← 0;
cy ← FALSE;
THROUGH [0..32) DO Step[]; ENDLOOP;
WHILE NOT BitOn[y.frac.highbits, HiBit] DO
-- normalize
Step[];
y.exp ← y.exp - 1;
ENDLOOP;
IF a.frac.lc # 0 THEN y.frac.lowbits ← Inline.BITOR[y.frac.lowbits, 1];
};
FRem: PUBLIC SAFE PROC [a, b: REAL, m: Mode] RETURNS [c: REAL] = TRUSTED { ERROR; };
END.
July 8, 1980 6:36 PM, L. Stewart; FRem dummy added
August 25, 1980 10:07 AM, L. Stewart; new Divide
January 28, 1981 11:50 AM, L. Stewart; Fix to Multiply (denormalized numbers)
August 27, 1982 11:24 am, L. Stewart, added TRUSTED everywhere