<> <> <> <> DIRECTORY Basics: TYPE USING [LongMult, LongNumber], PrincOpsUtils: TYPE USING [BITAND, BITOR, BITSHIFT, BITXOR], Real: FROM "IeeeFloat" USING [ Exception, ExceptionFlags, Extended, MinusInfinity, MinusZero, Mode, NoExceptions, NumberType, PlusInfinity, PlusZero, TrapNonTrappingNaN, TrapTrappingNaN]; IeeePack: PROGRAM IMPORTS Basics, PrincOpsUtils EXPORTS Real = { OPEN Basics, PrincOpsUtils, Real; RealException: PUBLIC SIGNAL [flags: ExceptionFlags] = CODE; <> Details: TYPE = MACHINE DEPENDENT RECORD [ sign: BOOL, sticky: BOOL, blank: [0..7777b], type: Real.NumberType]; Ext: TYPE = RECORD [ det: Details, exp: INTEGER, frac: LongNumber]; <> 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; <> 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]}; <> 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]}; <> 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]; <> 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; <> IF lo.lc>HalfLC OR (lo.lc=HalfLC AND BitOn[z.frac.lowbits, 1]) THEN { <> z.frac.lc _ z.frac.lc+1; IF z.frac.lc 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; <> <> 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]; <> RETURN[Pack[@y]]}; <> 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]; <> 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]}; <> 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}; <> 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 _ PrincOpsUtils.BITSHIFT[ PrincOpsUtils.BITAND[LN[r].highbits, ExponentMask], -ExponentShift]; z.exp _ z.exp - ExponentBias; z.det.type _ normal; z.frac.li _ LN[r].li; z.frac.highbits _ PrincOpsUtils.BITAND[HiFractionMask, z.frac.highbits]; SELECT z.exp FROM = DenormalizedExponent => { <> 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 => { <> 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 _ PrincOpsUtils.BITOR[HiddenBit, z.frac.highbits]}; }; <> Pack: PROC [z: POINTER TO Ext] RETURNS [r: REAL] = { trap: BOOL _ FALSE; i: Exception; IF thisTimeExceptions#NoExceptions THEN { <> 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]; <