
<<DragonRealIeeeA.mesa>>
<<Last Modified: Stewart January 14, 1984 4:49 pm>>
<<Last Modified: Paul Rovner May 4, 1983 9:53 am>>
<<Mesa implementation of floating point ops +, -, *, /, Rem>>

DIRECTORY
IeeeInternal,
PrincOpsUtils,
Real USING [
AddInfinityNaN, DivideInfinityNaN, MinusZero, MultiplyInfinityNaN,
NoExceptions, NumberType, PlusZero],
RealOps;

DragonRealIeeeA: PROGRAM IMPORTS IeeeInternal, PrincOpsUtils EXPORTS RealOps =
BEGIN OPEN IeeeInternal, PrincOpsUtils, 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 _ 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 _ 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
January 14, 1984 4:49 pm, Stewart, change to IeeeInternal