
    <<IeeeFloatA.mesa - Mesa implementation of floating point ops +, -, *, /, Rem>>
        <<Copyright © 1985 by Xerox Corporation.  All rights reserved.>>
        <<Stewart January 14, 1984 4:49 pm>>
        <<Rovner May 4, 1983 9:53 am>>
    <<Russ Atkinson (RRA) May 28, 1985 6:44:40 pm PDT>>

    DIRECTORY
        Ieee USING [ADC2, ADC3, BitOn, DeNormalize, Ext, funny, HiBit, HiddenBit, LongNumber, LowSignificandMask, Mul32, NormalType, Pack, 
        PostNormalize, RShift1in1, SetDivisionByZero, SetInvalidOperation, SignBit, StepTwo, TNnn, TNnz, TNzn, TNzz, Unpack],
        PrincOpsUtils USING [BITAND, BITOR, BITXOR],
        Real USING [AddInfinityNaN, DivideInfinityNaN, Fix, MinusZero, MultiplyInfinityNaN, NumberType, PlusZero],
        RealExceptions USING [ClearExceptions],
        RealOps USING [GetMode, Mode, SetMode];

    IeeeFloatA: PROGRAM
        IMPORTS Ieee, PrincOpsUtils, Real, RealExceptions, RealOps
        EXPORTS RealOps
        = BEGIN OPEN Ieee, PrincOpsUtils, Real, RealOps;

    NormalizedCardinal: PROC [g: CARDINAL] RETURNS [BOOL] = INLINE {
        RETURN [PrincOpsUtils.BITAND[g, Ieee.HiddenBit] # 0];
        };

    AddNormal: PROC [x, y: Ext] RETURNS [Ext, Ext] = {
        ediff: INTEGER _ x.exp - y.exp;
        normalized: BOOL;
        IF ediff > 0
            THEN DeNormalize[@y, -ediff]
            ELSE {DeNormalize[@x, ediff]; x.exp _ y.exp};
        normalized _ NormalizedCardinal[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 _ RShift1in1[x.frac];
                    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 _ RealOps.GetMode[].round = rm;
                        x.det.type _ zero;
                        }
                    ELSE {IF normalized THEN PostNormalize[@x]; };
                };
        IF y.det.sticky THEN x.det.sticky _ TRUE;
        RETURN [x, y];
        };

    FAdd: PUBLIC SAFE PROC [a, b: REAL, m: Mode] RETURNS [r: REAL] = TRUSTED {
        inner: PROC = {
            x: Ext _ Unpack[a];
            y: Ext  _ Unpack[b];
            {
                SELECT TRUE FROM
                    Basics.BITAND[funny, BITOR[LOOPHOLE[x.det], LOOPHOLE[y.det]]] = 0 => {
                        <<Normal case: Neither operand is infinity or NaN>>
                        SELECT NormalType[x.det, y.det] FROM
                            TNnn => NULL;
                            TNnz => GOTO RetA;
                            TNzn => GOTO RetB;
                            TNzz => {
                                a _ PlusZero;
                                IF m.round = rm AND x.det.sign AND y.det.sign THEN a _ MinusZero;
                                GO TO RetA;
                                };
                            ENDCASE => ERROR;
                        };
                    x.det.type = nan => GOTO RetX;
                    y.det.type = nan => GOTO RetY;
                    x.det.type = infinity AND y.det.type = infinity =>
                        IF m.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;
                [x, y] _ AddNormal[x, y];
                StepTwo[@x];
                GOTO RetX;
                EXITS
                    RetA => r _ a;
                    RetB => r _ b;
                    RetX => r _ Pack[@x];
                    RetY => r _ Pack[@y];
                };
            };
        DoWithMode[inner, m];
        };

    FSub: PUBLIC SAFE PROC [a, b: REAL, m: Mode] RETURNS [REAL] = TRUSTED {
        <<Negate b>>
        LOOPHOLE[b, LongNumber].highbits _ BITXOR[LOOPHOLE[b, LongNumber].highbits, SignBit];
        RETURN[FAdd[a, b, m]];
        };

    FMul: PUBLIC SAFE PROC [a, b: REAL, m: Mode] RETURNS [r: REAL] = TRUSTED {
        inner: PROC = {
            x: Ext _ Ieee.Unpack[a];
            y: Ext  _ Ieee.Unpack[b];
            {
                lo: LongNumber;
                x.det.sign _ x.det.sign # y.det.sign; -- due to different bias
                x.exp _ x.exp + y.exp + 1;
                SELECT TRUE FROM
                    Basics.BITAND[funny, BITOR[LOOPHOLE[x.det], LOOPHOLE[y.det]]] = 0 => {
                        <<Normal case: Neither operand is infinity or NaN>>
                        SELECT NormalType[x.det, y.det] FROM
                            TNnn => NULL;
                            TNnz => GOTO RetZero;
                            TNzn => GOTO RetZero;
                            TNzz => GOTO RetZero;
                            ENDCASE => ERROR;
                        };
                    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;
                [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 => r _ Pack[@x];
                    RetY => r _ Pack[@y];
                    RetZero => r _ IF x.det.sign THEN MinusZero ELSE PlusZero;
                };
            };
        DoWithMode[inner, m];
        };

    FDiv: PUBLIC SAFE PROC [a, b: REAL, m: Mode] RETURNS [c: REAL] = TRUSTED {
        inner: PROC = {
            x: Ext _ Unpack[a];
            y: Ext  _ Unpack[b];
            {
                x.det.sign _ x.det.sign # y.det.sign;
                {
                    SELECT TRUE FROM
                        Basics.BITAND[funny, BITOR[LOOPHOLE[x.det], LOOPHOLE[y.det]]] = 0 => {
                            <<Neither of the number types is infinity or NaN>>
                            SELECT NormalType[x.det, y.det] FROM
                                TNnn => {};
                                    <<this is usual flow>>
                                        <<RRA: this used to be the following bogus code:>>
                                            <<IF NormalizedCardinal[y.frac.highbits] THEN NULL ELSE GOTO ZeroDivide>>
                                TNnz => GOTO ZeroDivide;
                                TNzn => GOTO RetZero;
                                TNzz => GOTO Undef;
                                ENDCASE => ERROR;
                            };
                        x.det.type = nan => GOTO RetX;
                        y.det.type = nan => GOTO RetY;
                        x.det.type = infinity AND y.det.type = infinity => GO TO Undef;
                            <<both are some kind of infinity, so results are undefined>>
                        y.det.type = infinity => GOTO RetZero;
                        x.det.type = infinity => GOTO RetX;
                        ENDCASE => ERROR;
                            <<It should not be possible to get here!>>
                    EXITS
                        ZeroDivide => {
                            <<Divide by 0, so result is infnity>>
                            SetDivisionByZero[];
                            x.det.type _ infinity;
                            GOTO RetX;
                            };
                        Undef => {
                            <<The results are undefined (inf/inf or 0/0)>>
                            SetInvalidOperation[];
                            x.det.type _ nan;
                            x.frac.lc _ DivideInfinityNaN;
                            GOTO RetX;
                            };
                    };
                <<Both numbers need to be normalized for the Divide routine to work.>>
                IF NOT NormalizedCardinal[x.frac.highbits] THEN PostNormalize[@x];
                IF NOT NormalizedCardinal[y.frac.highbits] THEN PostNormalize[@y];
                x _ Divide[x, y];
                StepTwo[@x];
                GOTO RetX;
                EXITS
                    RetZero => c _ IF x.det.sign THEN MinusZero ELSE PlusZero;
                    RetX => c _ Pack[@x];
                    RetY => c _ Pack[@y];
                };
            };
        DoWithMode[inner, m];
        };

    Divide: PROC [a, b: Ext] RETURNS [y: Ext] = {
        cy: BOOL;
        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;
            };
        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 PrincOpsUtils.BITAND[y.frac.highbits, HiBit] = 0 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 {
        n: INT _ Real.Fix[a / b];
        c _ a - n * b;
        };

    DoWithMode: PROC [inner: PROC, mode: Mode] = {
        old: Mode _ RealOps.SetMode[mode];
        RealExceptions.ClearExceptions[];
        inner[ ! UNWIND => [] _ RealOps.SetMode[old]];
        [] _ RealOps.SetMode[old];
        };

    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 Ieee