! ---------------------------------------------------------------------------- ! FP: Core of floating point library ! ! Supplied for use with Inform 6 Serial number 020326 ! Release 1/2 ! (c) Kevin Bracey 2002 ! but freely usable (see manuals) ! ---------------------------------------------------------------------------- System_file; Constant FE_INVALID = $01; Constant FE_DIVBYZERO = $02; Constant FE_OVERFLOW = $04; Constant FE_UNDERFLOW = $08; Constant FE_INEXACT = $10; Constant FE_ALL_EXCEPT = $1F; Class FEnv with rounding FE_TONEAREST, status 0, printmode FE_PRINT_G, printprecision 6, inx_handler [; print "[** Floating-point trap: Inexact **]^"; quit; ], ufl_handler [; print "[** Floating-point trap: Underflow **]^"; quit; ], ofl_handler [; print "[** Floating-point trap: Overflow **]^"; quit; ], ivo_handler [; print "[** Floating-point trap: Invalid operation **]^"; quit; ], dvz_handler [; print "[** Floating-point trap: Division by zero **]^"; quit; ]; FEnv activefenv; ! Rounding modes Constant FE_TONEAREST = 1; Constant FE_UPWARD = 2; Constant FE_DOWNWARD = 3; Constant FE_TOWARDZERO = 4; ! Invalid operation reasons, stored in NaNs and passed to trap handlers Constant InvReason_SigNaN = 0; Constant InvReason_InitNaN = 1; Constant InvReason_MagSubInf = 4; Constant InvReason_InfTimes0 = 5; Constant InvReason_0TimesInf = 6; Constant InvReason_0Div0 = 7; Constant InvReason_InfDivInf = 8; Constant InvReason_InfRemX = 9; Constant InvReason_XRem0 = 10; Constant InvReason_SqrtNeg = 11; Constant InvReason_FixQNaN = 12; Constant InvReason_FixInf = 13; Constant InvReason_FixRange = 14; Constant InvReason_CompQNaN = 15; ! Destination format specifiers passed to trap handlers Constant FE_FMT_S = 0; ! Single Constant FE_FMT_X = 1; ! Extended Constant FE_FMT_P = 2; ! Packed decimal Constant FE_FMT_I = 3; ! Integer Constant FE_FMT_N = 4; ! None ! Operation codes passed to trap handlers Constant FE_OP_ADD = 0; Constant FE_OP_SUB = 1; Constant FE_OP_MUL = 2; Constant FE_OP_DIV = 3; Constant FE_OP_REM = 4; Constant FE_OP_CMP = 5; Constant FE_OP_CMPE = 6; Constant FE_OP_CONV = 7; Constant FE_OP_DEC = 8; Constant FE_OP_FIX = 9; Constant FE_OP_FLT = 10; Constant FE_OP_RND = 11; Constant FE_OP_SQRT = 12; Constant FE_OP_RAISE = 13; ! Decimal printing format Constant FE_PRINT_G = 0; Constant FE_PRINT_E = 1; Constant FE_PRINT_F = 2; ! Return values from fcmp[e][x] Constant FCMP_U = $8; Constant FCMP_L = $4; Constant FCMP_E = $2; Constant FCMP_G = $1; ! Global state variables Global fpstatus; ! Working copy of activefenv.fpstatus (for speed) ! Internal workings Array _fpscratch --> 6; Array _workF0 --> 3; Array _workF1 --> 3; Global _precision; Global _dest; Global _op; Global _rounding; Constant INXE = $1000; Constant UFLE = $0800; Constant OFLE = $0400; Constant DVZE = $0200; Constant IVOE = $0100; #Ifndef NULL; Constant NULL = $FFFF; #Endif; ! Internal utility functions [ _UCmp x y; @add x $8000 -> x; @add y $8000 -> y; @jg x y ?rtrue; @jl x y ?~rfalse; @ret -1; ]; [ _mul32 dest ah bh al bl m1 m2 tmp; ! Split a and b into two 8-bit halves @and ah $00FF -> al; @log_shift ah 0-8 -> ah; @and bh $00FF -> bl; @log_shift bh 0-8 -> bh; ! Four 8x8->16 multiplies m1 = al * bh; ! ( 0 m1h, m1l 0 ) m2 = ah * bl; ! + ( 0 m2h, m2l 0 ) ah = ah * bh; ! + ( ah , 0 ) al = al * bl; ! + ( 0 , al ) ! Add m1 into (ah,al) @log_shift m1 8 -> tmp; al = al + tmp; if (_UCmp(al,tmp)<0) ah++; @log_shift m1 0-8 -> tmp; ah = ah + tmp; ! Add m2 into (ah,al) @log_shift m2 8 -> tmp; al = al + tmp; if (_UCmp(al,tmp)<0) ah++; @log_shift m2 0-8 -> tmp; ah = ah + tmp; dest-->0 = ah; dest-->1 = al; ]; ! Clear specified status flags [ feclearexcept excepts; excepts = excepts & FE_ALL_EXCEPT; fpstatus = fpstatus &~ excepts; ]; ! Set specified status flags [ fesetexcept excepts; excepts = excepts & FE_ALL_EXCEPT; fpstatus = fpstatus | excepts; ]; ! Raise specified floating point exceptions [ feraiseexcept excepts round; round = fegetround(); if (excepts & FE_INVALID) { if (fpstatus & IVOE) activefenv.ivo_handler(NULL, FE_FMT_N, FE_OP_RAISE, round, InvReason_InitNaN); else fpstatus = fpstatus | FE_INVALID; } if (excepts & FE_DIVBYZERO) { if (fpstatus & DVZE) activefenv.dvz_handler(NULL, FE_FMT_N, FE_OP_RAISE, round); else fpstatus = fpstatus | FE_DIVBYZERO; } if (excepts & FE_OVERFLOW) { if (fpstatus & OFLE) activefenv.ofl_handler(NULL, FE_FMT_N, FE_OP_RAISE, round); else fpstatus = fpstatus | FE_OVERFLOW; } if (excepts & FE_UNDERFLOW) { if (fpstatus & UFLE) activefenv.ufl_handler(NULL, FE_FMT_N, FE_OP_RAISE, round); else fpstatus = fpstatus | FE_UNDERFLOW; } if (excepts & FE_INEXACT) { if (fpstatus & INXE) activefenv.inx_handler(NULL, FE_FMT_N, FE_OP_RAISE, round); else fpstatus = fpstatus | FE_INEXACT; } ]; ! Test specified status flags [ fetestexcept excepts; excepts = excepts & FE_ALL_EXCEPT; return fpstatus & excepts; ]; ! Get current rounding mode [ fegetround; return activefenv.rounding; ]; ! Set current rounding mode [ fesetround round; activefenv.rounding = round; ]; ! Store current floating-point environment in envp [ fegetenv envp; activefenv.fpstatus = fpstatus; FEnv.copy(envp, activefenv); ]; ! Store current floating-point environment in envp, clear all exceptions ! and disable traps [ feholdexcept envp; fegetenv(envp); fpstatus = 0; ]; ! Restore floating-point environment from envp [ fesetenv envp; Fenv.copy(activefenv, envp); fpstatus = activefenv.fpstatus; ]; ! Remember exception status, restore environment, then raise the saved ! exceptions [ feupdateenv envp oldstatus; oldstatus = fetestexcept(FE_ALL_EXCEPT); fesetenv(envp); feraiseexcept(oldstatus); ]; ! Test specified trap enables [ fetesttrap traps; traps = traps & FE_ALL_EXCEPT; @log_shift fpstatus 0-8 -> sp; @and sp traps -> sp; @ret_popped; ]; ! Enable specified traps [ feenabletrap traps; traps = traps & FE_ALL_EXCEPT; @log_shift traps 8 -> sp; @or fpstatus sp -> fpstatus; ]; ! Disable specified traps [ fedisabletrap traps; traps = traps & FE_ALL_EXCEPT; @log_shift traps 8 -> traps; fpstatus = fpstatus &~ traps; ]; ! Set the trap handlers [ fesettraphandler routine traps; if (traps & FE_INVALID) activefenv.ivo_handler = routine; if (traps & FE_DIVBYZERO) activefenv.dvz_handler = routine; if (traps & FE_OVERFLOW) activefenv.ofl_handler = routine; if (traps & FE_UNDERFLOW) activefenv.ufl_handler = routine; if (traps & FE_INEXACT) activefenv.inx_handler = routine; ]; [ fegettraphandler trap; switch (trap) { FE_INVALID: return activefenv.ivo_handler; FE_DIVBYZERO: return activefenv.dvz_handler; FE_OVERFLOW: return activefenv.ofl_handler; FE_UNDERFLOW: return activefenv.ufl_handler; FE_INEXACT: return activefenv.inx_handler; default: return NULL; } ]; [ fesetprintmode mode oldmode; oldmode = activefenv.printmode; activefenv.printmode = mode; return oldmode; ]; [ fesetprintprecision prec oldprec; oldprec = activefenv.printprecision; activefenv.printprecision = prec; return oldprec; ]; [ isnan x h; h = x-->0; @test h $7F80 ?~rfalse; @and h $007F -> sp; @loadw x 1 -> sp; @or sp sp -> sp; @jz sp ?rfalse; rtrue; ]; [ isnanx x; @loadw x 0 -> sp; @test sp $07FF ?~rfalse; @loadw x 1 -> sp; @loadw x 2 -> sp; @or sp sp -> sp; @jz sp ?rfalse; rtrue; ]; [ _isnan sex mhi mlo; @test sex $07FF ?~rfalse; @or mhi mlo -> sp; @jz sp ?rfalse; rtrue; ]; [ issignalling x h; @loadw x 0 -> h; @and h $7FC0 -> sp; @je sp $7F80 ?~rfalse; @and h $007F -> sp; @loadw x 1 -> sp; @or sp sp -> sp; @jz sp ?rfalse; rtrue; ]; [ issignallingx x h; @loadw x 0 -> sp; @test sp $07FF ?~rfalse; @loadw x 1 -> h; @loadw x 2 -> sp; @or h sp -> sp; @jz sp ?rfalse; @test h $4000 ?rfalse; rtrue; ]; [ _issignalling sex mhi mlo; @test sex $07FF ?~rfalse; @or mhi mlo -> sp; @jz sp ?rfalse; @test mhi $4000 ?rfalse; rtrue; ]; [ isinf x; @loadw x 1 -> sp; @jz sp ?~rfalse; @loadw x 0 -> sp; @and sp $7FFF -> sp; @je sp $7F80 ?~rfalse; rtrue; ]; [ isinfx x; @loadw x 1 -> sp; @jz sp ?~rfalse; @loadw x 2 -> sp; @jz sp ?~rfalse; @loadw x 0 -> sp; @test sp $07FF ?~rfalse; rtrue; ]; [ isfinite x; @loadw x 0 -> sp; @test x $7F80 ?~rtrue; rfalse; ]; [ isfinitex x; @loadw x 0 -> sp; @test x $07FF ?~rtrue; rfalse; ]; ! Return values from fpclassify[x] Constant FP_NANS = $0001; Constant FP_NANQ = $0002; Constant FP_ZEROM = $0004; Constant FP_ZEROP = $0008; Constant FP_SUBNORMALM = $0010; Constant FP_SUBNORMALP = $0020; Constant FP_NORMALM = $0040; Constant FP_NORMALP = $0080; Constant FP_INFINITYM = $0100; Constant FP_INFINITYP = $0200; Constant FP_NAN = FP_NANS | FP_NANQ; Constant FP_ZERO = FP_ZEROM | FP_ZEROP; Constant FP_SUBNORMAL = FP_SUBNORMALM | FP_SUBNORMALP; Constant FP_NORMAL = FP_NORMALM | FP_NORMALP; Constant FP_INFINITY = FP_INFINITYM | FP_INFINITYP; [ fpclassify x; return fpclassifyx(_Promote(x)); ]; [ fpclassifyx x s h l e; s = x-->0; h = x-->1; l = x-->2; e = s & $07FF; if (e == $07FF) { if ((h|l) == 0) x = FP_INFINITYM; else if (h & $4000) return FP_NANQ; else return FP_NANS; } else if ((h|l) == 0) x = FP_ZEROM; else if (h >= 0) x = FP_SUBNORMALM; else x = FP_NORMALM; if (s >= 0) @log_shift x 1 -> x; return x; ]; [ fcmp_common OP1 OP2 OP1emh OP2emh OP1sxm OP2sxm OP1mlo OP2mlo; OP1sxm = OP1-->0; OP1mlo = OP1-->1; OP2sxm = OP2-->0 + $8000; OP2mlo = OP2-->1; OP1emh = OP1sxm & $7FFF; OP2emh = OP2sxm & $7FFF; @jg OP1emh OP2emh ?op1bigger; @jl OP1emh OP2emh ?op2bigger; OP1 = _UCmp(OP1mlo, OP2mlo); @jg OP1 0 ?op1bigger; @jz OP1 ?equalmag; .op2bigger; if (OP2sxm>=0) return FCMP_G; else return FCMP_L; .op1bigger; if (OP1sxm>=0) return FCMP_G; else return FCMP_L; .equalmag; if ((OP1sxm<0 && OP2sxm>=0) || (OP1sxm>=0 && OP2sxm<0) || (OP1emh|OP1mlo)==0) return FCMP_E; if (OP1sxm>=0) return FCMP_G; else return FCMP_L; ]; [ fcmpe OP1 OP2; if (isnan(OP1) || isnan(OP2)) return _RaiseCmpIVO(FE_FMT_S, FE_OP_CMPE, OP1, OP2); return fcmp_common(OP1,OP2); ]; [ fcmp OP1 OP2; if (isnan(OP1) || isnan(OP2)) { if (issignalling(OP1) || issignalling(OP2)) return _RaiseCmpIVO(FE_FMT_S, FE_OP_CMP, OP1, OP2); else return FCMP_U; } return fcmp_common(OP1,OP2); ]; [ fcmpx_common OP1 OP2 OP1exp OP2exp OP1sex OP2sex OP1mhi OP2mhi OP1mlo OP2mlo; OP1sex = OP1-->0; OP1mhi = OP1-->1; OP1mlo = OP1-->2; OP2sex = OP2-->0 + $8000; OP2mhi = OP2-->1; OP2mlo = OP2-->2; OP1exp = OP1sex & $7FFF; OP2exp = OP2sex & $7FFF; @jg OP1exp OP2exp ?op1bigger; @jl OP1exp OP2exp ?op2bigger; OP1 = _UCmp(OP1mhi, OP2mhi); @jg OP1 0 ?op1bigger; @jl OP1 0 ?op2bigger; OP1 = _UCmp(OP1mlo, OP2mlo); @jg OP1 0 ?op1bigger; @jz OP1 ?equalmag; .op2bigger; if (OP2sex>=0) return FCMP_G; else return FCMP_L; .op1bigger; if (OP1sex>=0) return FCMP_G; else return FCMP_L; .equalmag; if ((OP1sex<0 && OP2sex>=0) || (OP1sex>=0 && OP2sex<0) || (OP1exp|OP1mhi|OP1mlo)==0) return FCMP_E; if (OP1sex>=0) return FCMP_G; else return FCMP_L; ]; [ fcmpex OP1 OP2; if (isnanx(OP1) || isnanx(OP2)) return _RaiseCmpIVO(FE_FMT_X, FE_OP_CMPE, OP1, OP2); return fcmpx_common(OP1,OP2); ]; [ fcmpx OP1 OP2; if (isnanx(OP1) || isnanx(OP2)) { if (issignallingx(OP1) || issignallingx(OP2)) return _RaiseCmpIVO(FE_FMT_X, FE_OP_CMP, OP1, OP2); else return FCMP_U; } return fcmpx_common(OP1,OP2); ]; [ feq OP1 OP2; return fcmp (OP1,OP2) == FCMP_E; ]; [ fne OP1 OP2; return fcmp (OP1,OP2) ~= FCMP_E; ]; [ fgt OP1 OP2; return fcmpe (OP1,OP2) == FCMP_G; ]; [ fge OP1 OP2; return fcmpe (OP1,OP2) & (FCMP_G|FCMP_E); ]; [ flt OP1 OP2; return fcmpe (OP1,OP2) == FCMP_L; ]; [ fle OP1 OP2; return fcmpe (OP1,OP2) & (FCMP_L|FCMP_E); ]; [ funordered OP1 OP2; return fcmp (OP1,OP2) == FCMP_U; ]; [ flg OP1 OP2; return fcmpe (OP1,OP2) & (FCMP_L|FCMP_G); ]; [ fleg OP1 OP2; return fcmpe (OP1,OP2) ~= FCMP_U; ]; [ fug OP1 OP2; return fcmp (OP1,OP2) & (FCMP_U|FCMP_G); ]; [ fuge OP1 OP2; return fcmp (OP1,OP2) ~= FCMP_L; ]; [ ful OP1 OP2; return fcmp (OP1,OP2) & (FCMP_U|FCMP_L); ]; [ fule OP1 OP2; return fcmp (OP1,OP2) ~= FCMP_G; ]; [ fue OP1 OP2; return fcmp (OP1,OP2) & (FCMP_U|FCMP_E); ]; [ feqx OP1 OP2; return fcmpx (OP1,OP2) == FCMP_E; ]; [ fnex OP1 OP2; return fcmpx (OP1,OP2) ~= FCMP_E; ]; [ fgtx OP1 OP2; return fcmpex(OP1,OP2) == FCMP_G; ]; [ fgex OP1 OP2; return fcmpex(OP1,OP2) & (FCMP_G|FCMP_E); ]; [ fltx OP1 OP2; return fcmpex(OP1,OP2) == FCMP_L; ]; [ flex OP1 OP2; return fcmpex(OP1,OP2) & (FCMP_L|FCMP_E); ]; [ funorderedx OP1 OP2; return fcmpx (OP1,OP2) == FCMP_U; ]; [ flgx OP1 OP2; return fcmpex(OP1,OP2) & (FCMP_L|FCMP_G); ]; [ flegx OP1 OP2; return fcmpex(OP1,OP2) ~= FCMP_U; ]; [ fugx OP1 OP2; return fcmpx (OP1,OP2) & (FCMP_U|FCMP_G); ]; [ fugex OP1 OP2; return fcmpx (OP1,OP2) ~= FCMP_L; ]; [ fulx OP1 OP2; return fcmpx (OP1,OP2) & (FCMP_U|FCMP_L); ]; [ fulex OP1 OP2; return fcmpx (OP1,OP2) ~= FCMP_G; ]; [ fuex OP1 OP2; return fcmpx (OP1,OP2) & (FCMP_U|FCMP_E); ]; [ fmax dest OP1 OP2 cmp; cmp = fcmp(OP1,OP2); if (cmp & (FCMP_E|FCMP_G)) jump ret1; if (cmp & FCMP_L) jump ret2; if (isnan(OP1)) jump ret2; .ret1; @copy_table OP1 dest 4; return; .ret2; @copy_table OP2 dest 4; return; ]; [ fmaxx dest OP1 OP2 cmp; cmp = fcmpx(OP1,OP2); if (cmp & (FCMP_E|FCMP_G)) jump ret1; if (cmp & FCMP_L) jump ret2; if (isnanx(OP1)) jump ret2; .ret1; @copy_table OP1 dest 6; return; .ret2; @copy_table OP2 dest 6; return; ]; [ fmin dest OP1 OP2 cmp; cmp = fcmp(OP1,OP2); if (cmp & (FCMP_E|FCMP_L)) jump ret1; if (cmp & FCMP_G) jump ret2; if (isnan(OP1)) jump ret2; .ret1; @copy_table OP1 dest 4; return; .ret2; @copy_table OP2 dest 4; return; ]; [ fminx dest OP1 OP2 cmp; cmp = fcmpx(OP1,OP2); if (cmp & (FCMP_E|FCMP_L)) jump ret1; if (cmp & FCMP_G) jump ret2; if (isnanx(OP1)) jump ret2; .ret1; @copy_table OP1 dest 6; return; .ret2; @copy_table OP2 dest 6; return; ]; ! No need for rounding as all FP formats ! can hold a 16-bit int. [ fitos dst n; _float(FE_FMT_S, dst, n); ]; [ fitox dst n; _float(FE_FMT_X, dst, n); ]; [ _float prec dst n sign exp; _op = FE_OP_FLT; _precision = prec; _rounding = 0; _dest = dst; if (n < 0) { sign = $8000; n = -n; } if (n == 0) { exp = 0; } else { n = _Normalise(1023+15, n); exp = n-->0; n = n-->1; } _RoundResult(sign, exp, n); ]; [ fstoi src rnd; return fxtoi(_Promote(src), rnd); ]; [ fxtoi tmp rnd OPsex OPmhi OPmlo exp mhi mlo res dir; _rounding = rnd; !print "fxtoi(", (fhexx) tmp, ")^"; OPsex = tmp-->0; OPmhi = tmp-->1; OPmlo = tmp-->2; exp = OPsex & $07FF; ! Want to slide the binary point to the bottom tmp = 1023 + 31 - exp; if (tmp <= 0) jump outofrange; _precision = FE_FMT_X; res = _Denorm(OPmhi, OPmlo, tmp); mhi = res-->0; mlo = res-->1; tmp = res-->2; res = _RoundNum(OPsex & $8000, exp, mhi, mlo, tmp); mhi = res-->0; mlo = res-->1; dir = res-->4; ! Now have a 32-bit number in mhi, mlo if (OPsex < 0) { ! 2's complement it mhi = ~mhi; mlo = -mlo; if (mlo == 0) ++mhi; } if ((mlo >= 0 && mhi ~= 0) || (mlo < 0 && mhi ~= -1)) jump outofrange; if (dir) { if (fpstatus & INXE) mlo = activefenv.inx_handler(mlo, FE_FMT_I, FE_OP_FIX, _rounding); else fpstatus = fpstatus | FE_INEXACT; } return mlo; .outofrange; return _RaiseFixIVO(OPsex, OPmhi, OPmlo); ]; [ fstox dst src noexcept sign exp mhi mlo res; mhi = src-->0; mlo = src-->1; ! print "fstox(", (hex) mhi, (hex) mlo, ")="; exp = (mhi & $7F80); @log_shift exp 0-7 -> exp; sign = mhi & $8000; mhi = mhi & $007F; @log_shift mhi 8 -> mhi; @log_shift mlo 0-8 -> sp; @or mhi sp -> mhi; @log_shift mlo 8 -> mlo; if (exp == 0) { if (mhi | mlo) { ! Subnormal res = _Normalise(1023-126, mhi, mlo); exp = res-->0; mhi = res-->1; mlo = res-->2; } else { ! Zero } } else if (exp == $FF) { ! Infinite or NaN exp = $7FF; if (mhi | mlo) { if ((~~noexcept) && (mhi & $4000)==0) { ! Conversion of signalling NaN _dest = dst; _precision = FE_FMT_X; _op = FE_OP_CONV; _RaiseIVO(InvReason_SigNaN, sign | exp, mhi, mlo); return; } } } else { ! Normal mhi = mhi | $8000; exp = exp + (1023 - 127); } dst-->0 = sign | exp; dst-->1 = mhi; dst-->2 = mlo; ! print (hex) dst-->0, "|", (hex) dst-->1, (hex) dst-->2, "^"; ]; [ fxtos dst src rnd sign exp sex mhi mlo; sex = src-->0; mhi = src-->1; mlo = src-->2; !print "fstox(", (hex) mhi, (hex) mlo, ")="; exp = sex & $07FF; sign = sex & $8000; _dest = dst; _precision = FE_FMT_S; _rounding = rnd; _op = FE_OP_CONV; _RoundResult(sign, exp, mhi, mlo); ]; [ fcpy dst src; @copy_table src dst 4; ]; [ fcpyx dst src; @copy_table src dst 6; ]; [ fneg dst src; dst-->1 = src-->1; dst-->0 = src-->0 + $8000; ]; [ fnegx dst src; dst-->2 = src-->2; dst-->1 = src-->1; dst-->0 = src-->0 + $8000; ]; [ fabs dst src; dst-->1 = src-->1; dst-->0 = src-->0 & $7FFF; ]; [ fabsx dst src; dst-->2 = src-->2; dst-->1 = src-->1; dst-->0 = src-->0 & $7FFF; ]; [ _Denorm h l s w b t1 grs; !print "Denormalising ", (hex) h, (hex) l, " by ", s, " bits^"; @log_shift s 0-4 -> w; ! words to shift b = s & $F; ! bits to shift (0-15) ! Can't guarantee that x << 16 == 0, so must skirt around that case if (b ~= 0) { t1 = 16 - b; b = -b; @log_shift l t1 -> grs; ! bottom b bits of l into grs @log_shift l b -> l; ! shift l down b bits @log_shift h t1 -> t1; ! bottom b bits of h into l l = l | t1; @log_shift h b -> h; ! shift h down b bits } if (w == 1 || w >= 3) { if (grs) grs = 1; grs = l | grs; l = h; h = 0; } if (w >= 2) { if (grs || l) grs = 1; grs = h | grs; l = 0; h = 0; } ! print "Result is ", (hex) h, (hex) l, "/", (hex) grs; new_line; _fpscratch-->0=h; _fpscratch-->1=l; _fpscratch-->2=grs; return _fpscratch; ]; ! Normalise extended precision number - deals with weird zeros, unnormalised ! operands and infinities properly. Useful for coercing the result of _Promote ! back into something suitable for user (eg trap handler) consumption. ! Assumes the parameter is representable in standard extended format without ! rounding (eg was originally a user-supplied standard or extended number). [ _fnrmx dest OP sex mhi mlo exp; sex = OP-->0 & $87FF; ! strip out nasty bits (just in case) mhi = OP-->1; mlo = OP-->2; exp = sex & $07FF; @copy_table OP dest 6; if (exp == $7FF) ! NaNs and infinities OK { dest-->1 = mhi & $7FFF; ! but ensure J clear return; } if ((mhi|mlo)==0) ! Check for zeros { dest-->0 = sex & $8000; ! Ensure exponent is 0 return; } if (mhi < 0) ! Already normalised return; ! Normalise it OP = _Normalise(exp, mhi, mlo); exp = OP-->0; mhi = OP-->1; mlo = OP-->2; ! If subnormal, denormalise again if (exp < 0) { OP = _Denorm(mhi, mlo, -exp); mhi = OP-->0; mlo = OP-->1; exp = 0; } dest-->0 = (sex & $8000) | exp; dest-->1 = mhi; dest-->2 = mlo; ]; ! Normalise (mhi,mlo), by shifting bits up such that the MSB of mhi is set. ! Returns adjusted (possibly negative) exponent. (mhi,mlo) may be zero or ! already normal (in which case no change is made). [ _Normalise exp mhi mlo tmp; ! print "Normalising ", (hex) mhi, (hex) mlo, (char) '/', exp; new_line; if (mhi == 0) { if (mlo == 0) jump out; mhi = mlo; mlo = 0; exp = exp - 16; } if ((mhi & $FF00) == 0) { @log_shift mhi 8 -> mhi; tmp = 8; } if ((mhi & $F000) == 0) { @log_shift mhi 4 -> mhi; tmp = tmp + 4; } if ((mhi & $C000) == 0) { @log_shift mhi 2 -> mhi; tmp = tmp + 2; } if (mhi >= 0) { mhi = mhi + mhi; tmp ++; } @sub tmp 16 -> sp; @log_shift mlo sp -> sp; @or mhi sp -> mhi; @log_shift mlo tmp -> mlo; exp = exp - tmp; .out; ! print "Result is ", (hex) mhi, (hex) mlo, (char) '/', exp; new_line; _fpscratch-->0 = exp; _fpscratch-->1 = mhi; _fpscratch-->2 = mlo; return _fpscratch; ]; [ _RoundNum RNDsgn RNDexp RNDmhi RNDmlo RNDgrs RNDdir dir lsb; if (_precision == FE_FMT_S) { if (RNDgrs) RNDgrs = 1; @log_shift RNDmlo 8 -> sp; @or sp RNDgrs -> RNDgrs; RNDmlo = RNDmlo & $FF00; lsb = $0100; } else lsb = $0001; if (_rounding == 0) _rounding = activefenv.rounding; switch (_rounding) { FE_TONEAREST: if (RNDgrs < 0) ! round (top) bit set { if (RNDgrs ~= $8000) ! sticky bits set dir = 1; else ! halfway case { if (RNDdir < 0 || (RNDdir == 0 && (RNDmlo & lsb))) dir = 1; else dir = -1; } } else if (RNDgrs > 0) dir = -1; FE_DOWNWARD: if (RNDgrs) dir = -(RNDsgn+1); ! = +32767 or -1 FE_UPWARD: if (RNDgrs) dir = RNDsgn + 1; ! = -32767 or +1 FE_TOWARDZERO: if (RNDgrs) dir = -1; } if (dir > 0) { RNDmlo = RNDmlo + lsb; if (RNDmlo == 0) { if (++RNDmhi == $0) ! Mantissa overflow { RNDmhi = $8000; RNDexp++; } } } ! Update rounding so far if (dir) RNDdir = dir; _fpscratch-->0 = RNDmhi; _fpscratch-->1 = RNDmlo; _fpscratch-->2 = RNDexp; _fpscratch-->3 = dir; ! direction of this rounding _fpscratch-->4 = RNDdir; return _fpscratch; ]; [ _ReturnResult sex mhi mlo tmp tmp2; if (_precision == FE_FMT_X) { _dest-->0 = sex; _dest-->1 = mhi; _dest-->2 = mlo; } else { ! Simple narrowing - input should be either: ! a) finite with exponent biased for single, unnormalised if necessary ! b) infinity/NaN with exponent = $7FFF tmp = sex & $FF; @log_shift tmp 7 -> tmp; mhi = mhi & $7FFF; @log_shift mhi 0-8 -> tmp2; _dest-->0 = (sex & $8000) | tmp | tmp2; @log_shift mhi 8 -> mhi; @log_shift mlo 0-8 -> mlo; _dest-->1 = mhi | mlo; } ]; ! "Exact", normalised result provided, as extended number split into 5 parts. ! Round it to destination precision, then check for over/underflow. ! Denormalise if necessary, and store. [ _RoundResult RNDsgn RNDexp RNDmhi RNDmlo RNDgrs RNDdir ExpMin ExpMax BiasAdjust res; !print "_RoundResult(",RNDsgn, (hex)RNDexp, " "; !print (hex) RNDmhi, (hex) RNDmlo, "|", (hex) RNDgrs, (char) ')'; new_line; res=_RoundNum(RNDsgn, RNDexp, RNDmhi, RNDmlo, RNDgrs, RNDdir); RNDmhi = res-->0; RNDmlo = res-->1; RNDexp = res-->2; RNDdir = res-->4; ! Rebias exponent to destination format if (_precision == FE_FMT_S) { RNDexp = RNDexp + 127 - 1023; ExpMin = $01; ExpMax = $FE; BiasAdjust = 192; } else { ExpMin = $000; ExpMax = $7FE; BiasAdjust = 1536; } if (RNDexp < ExpMin || RNDexp > ExpMax) { ! print "Exponent out of range^"; if (RNDexp < ExpMin) { ! Potential underflow if (RNDmhi | RNDmlo) { if (fpstatus & UFLE) { ! Take underflow trap if (RNDexp + BiasAdjust < ExpMin) { ! Massive underflow if (_precision == FE_FMT_S) BiasAdjust = 127 - RNDexp; else BiasAdjust = 1023 - RNDexp; } BiasAdjust = -BiasAdjust; res = activefenv.ufl_handler; jump oflufl; } else { res = _Denorm(RNDmhi, Rndmlo, ExpMin - RNDexp); RNDmhi = res-->0; RNDmlo = res-->1; RNDgrs = res-->2; RNDexp = 0; res = _RoundNum(RNDsgn, RNDexp, RNDmhi, RNDmlo, RNDgrs, RNDdir); RNDmhi = res-->0; RNDmlo = res-->1; RNDdir = res-->4; ! Check it didn't round back up to be normalised if (RNDmhi < 0) RNDexp = ExpMin; if (res-->3) ! Denormalisation loss { fpstatus = fpstatus | FE_UNDERFLOW; ! print "Underflow (denormalisation loss)^"; } } } else RNDexp = 0; } else ! RNDexp > ExpMax { if (fpstatus & OFLE) { if (RNDexp - BiasAdjust > ExpMax) { ! Massive overflow if (_precision == FE_FMT_S) BiasAdjust = RNDexp - 127; else BiasAdjust = RNDexp - 1023; } res = activefenv.ofl_handler; .oflufl; _ReturnResult(RNDsgn | (RNDexp - BiasAdjust), RNDmhi, RNDmlo); @call_vn2 res _dest _precision _op _rounding BiasAdjust; return; } else { fpstatus = fpstatus | FE_OVERFLOW; res = _rounding; if (res == 0) res = fegetround(); if (res == FE_TONEAREST || (res == FE_UPWARD && RNDsgn >= 0) || (res == FE_DOWNWARD && RNDsgn < 0)) { RNDexp = $7FF; RNDmhi = 0; RNDmlo = 0; RNDdir = 1; } else { RNDexp = ExpMax; RNDmhi = $FFFF; RNDmlo = $FFFF; RNDdir = -1; } } } } !print (hex) RNDexp, (hex) RNDmhi, (hex) RNDmlo; new_line; _ReturnResult(RNDsgn | RNDexp, RNDmhi, RNDmlo); if (RNDdir ~= 0) { if (fpstatus & INXE) activefenv.inx_handler(_dest, _precision, _op, _rounding); else fpstatus = fpstatus | FE_INEXACT; } ]; ! Fast, non-excepting promotion of a number from single to extended ! for internal use. Subnormal numbers will be left unnormalised, ! zeros will have unusual exponents. [ _Promote OP sex mhi mlo exp; sex = OP-->0; mlo = OP-->1; exp = sex & $7F80; mhi = sex & $007F; @log_shift exp 0-7 -> exp; @log_shift mhi 8 -> mhi; @log_shift mlo 0-8 -> sp; @or mhi sp -> mhi; @log_shift mlo 8 -> mlo; if (exp == 0) exp = 1023-126; else if (exp == $FF) exp = $7FF; else { exp = exp + (1023-127); mhi = mhi | $8000; } _fpscratch-->0 = (sex & $8000) | exp; _fpscratch-->1 = mhi; _fpscratch-->2 = mlo; return _fpscratch; ]; [ _RaiseFixIVO sex mhi mlo reason; if (fpstatus & IVOE) { if ((sex & $07FF) == $07FF) { if (mhi|mlo) { if (mhi & $4000) reason = InvReason_FixQNaN; else reason = InvReason_SigNaN; } else reason = InvReason_FixInf; } else reason = InvReason_FixRange; _workF0-->0 = sex; _workF0-->1 = mhi; _workF0-->2 = mlo; _fnrmx(_workF0, _workF0); print (frawx) _workF0; new_line; sex = activefenv.ivo_handler; if (_rounding == 0) _rounding = fegetround(); @call_vs2 sex 0 FE_FMT_I FE_OP_FIX _rounding reason _workF0 -> sp; @ret_popped; } else { fpstatus = fpstatus | FE_INVALID; ! We return maximum or minimum integer, depending on sign if (sex >= 0) return $7FFF; else return $8000; } ]; [ _RaiseCmpIVO fmt op OP1 OP2 reason; if (fpstatus & IVOE) { if (fmt == FE_FMT_S) OP1 = _Promote(OP1); _fnrmx(_workF0, OP1); if (fmt == FE_FMT_S) OP2 = _Promote(OP2); _fnrmx(_workF1, OP2); if (((OP1-->1 & OP2-->1) & $4000) == 0) reason = InvReason_SigNaN; else reason = InvReason_CompQNaN; fmt = activefenv.ivo_handler; @call_vs2 fmt 0 FE_FMT_I op 0 reason _workF0 _workF1 -> sp; @ret_popped; } else { fpstatus = fpstatus | FE_INVALID; return FCMP_U; } ]; [ _RaiseIVO reason OP1sex OP1mhi OP1mlo OP2sex OP2mhi OP2mlo tmp; if (fpstatus & IVOE) { tmp = activefenv.ivo_handler; _workF0-->0 = OP1sex; _workF0-->1 = OP1mhi; _workF0-->2 = OP1mlo; _fnrmx(_workF0, _workF0); @check_arg_count 7 ?haveop2; @call_vn2 tmp _dest _precision _op _rounding reason _workF0; return; .haveop2; _workF1-->0 = OP2sex; _workF1-->1 = OP2mhi; _workF1-->2 = OP2mlo; _fnrmx(_workF1, _workF1); if (_rounding == 0) _rounding = fegetround(); @call_vn2 tmp _dest _precision _op _rounding reason _workF0 _workF1; } else { fpstatus = fpstatus | FE_INVALID; @log_shift reason 8 -> reason; _ReturnResult($07FF, $4000, reason); } ]; [ _RaiseDVZ OP1sex OP1mhi OP1mlo OP2sex tmp; if (fpstatus & DVZE) { _workF0-->0 = OP1sex; _workF0-->1 = OP1mhi; _workF0-->2 = OP1mlo; _fnrmx(_workF0, _workF0); _workF1-->0 = OP2sex & $8000; _workF1-->1 = $0000; _workF1-->2 = $0000; tmp = activefenv.dvz_handler; if (_rounding == 0) _rounding = fegetround(); @call_vn2 tmp _dest _precision _op _rounding 0 _workF0 _workF1; } else { fpstatus = fpstatus | FE_DIVBYZERO; ! Return appropriately signed infinity _ReturnResult($07FF | ((OP1sex & $8000) + (OP2sex & $8000))); } ]; [ _dyadic op func prec dest OP1 OP2 rnd OP1sex OP1mhi OP1mlo OP2sex OP2mhi OP2mlo; _op = op; _precision = prec; _dest = dest; _rounding = rnd; if (prec==0) OP1 = _Promote(OP1); OP1sex = OP1-->0; OP1mhi = OP1-->1; OP1mlo = OP1-->2; if (prec==0) OP2 = _Promote(OP2); OP2sex = OP2-->0; OP2mhi = OP2-->1; OP2mlo = OP2-->2; @test OP1sex $07FF ?uncommon; @test OP2sex $07FF ?uncommon; .proceed; @call_vn2 func OP1sex OP1mhi OP1mlo OP2sex OP2mhi OP2mlo; return; .uncommon; if (_issignalling(OP1sex, OP1mhi, OP1mlo) || _issignalling(OP2sex, OP2mhi, OP2mlo)) { _RaiseIVO(InvReason_SigNaN, OP1sex, OP1mhi, OP1mlo, OP2sex, OP2mhi, OP2mlo); return; } else if (_isnan(OP1sex, OP1mhi, OP1mlo)) { _ReturnResult(OP1sex, OP1mhi, OP1mlo); return; } else if (_isnan(OP2sex, OP2mhi, OP2mlo)) { _ReturnResult(OP2sex, OP2mhi, OP2mlo); return; } jump proceed; ]; [ fadd dest OP1 OP2 rnd; return _dyadic(FE_OP_ADD, _fadd, FE_FMT_S, dest, OP1, OP2, rnd); ]; [ faddx dest OP1 OP2 rnd; _dyadic(FE_OP_ADD, _fadd, FE_FMT_X, dest, OP1, OP2, rnd); ]; [ fsub dest OP1 OP2 rnd; _dyadic(FE_OP_SUB, _fsub, FE_FMT_S, dest, OP1, OP2, rnd); ]; [ fsubx dest OP1 OP2 rnd; _dyadic(FE_OP_SUB, _fsub, FE_FMT_X, dest, OP1, OP2, rnd); ]; [ _fsub OP1sex OP1mhi OP1mlo OP2sex OP2mhi OP2mlo; _fadd(OP1sex, OP1mhi, OP1mlo, OP2sex+$8000, OP2mhi, OP2mlo); ]; [ _fadd OP1sex OP1mhi OP1mlo OP2sex OP2mhi OP2mlo RNDgrs RNDexp tmp tmp2 OP1exp OP2exp; OP1exp = OP1sex & $07FF; OP2exp = OP2sex & $07FF; if (OP1exp == $7FF || OP2exp == $7FF) { if (OP2exp ~= $7FF) ! infinity + finite { _ReturnResult(OP1sex, OP1mhi, OP1mlo); return; } if (OP1exp ~= $7FF) ! finite + infinity { _ReturnResult(OP2sex, OP2mhi, OP2mlo); return; } ! two infinities if ((OP1sex & $8000) == (OP2sex & $8000)) { _ReturnResult(OP1sex, OP1mhi, OP1mlo); return; } else { if (_op == FE_OP_SUB) OP2sex = OP2sex + $8000; _RaiseIVO(InvReason_MagSubInf, OP1sex, OP1mhi, OP1mlo, OP2sex, OP2mhi, OP2mlo); return; } } ! We let zeros and subnormals through. Algorithm will work ! fine, but we may end up with a subnormal result. ! Denormalise the smaller operand to the same exponent ! as the larger. if (OP1exp < OP2exp) { RNDexp = OP2exp; tmp = _Denorm(OP1mhi, OP1mlo, OP2exp - OP1exp); OP1mhi = tmp-->0; OP1mlo = tmp-->1; tmp = tmp-->2; tmp2 = 0; } else if (OP1exp > OP2exp) { RNDexp = OP1exp; tmp = _Denorm(OP2mhi, OP2mlo, OP1exp - OP2exp); OP2mhi = tmp-->0; OP2mlo = tmp-->1; tmp2 = tmp-->2; tmp = 0; } else { RNDexp = OP1exp; tmp = tmp2 = 0; } ! Don't need original numbers any longer OP1sex = OP1sex & $8000; OP2sex = OP2sex & $8000; ! Now OP1sex/OP2sex = signs of OP1/OP2 ! OP1mhi/OP1mlo = operand 1 mantissa ! RNDexp = prospective result exponent ! OP2mhi/OP2mlo = operand 2 mantissa ! tmp = operand 1 guard, round and sticky bits ! tmp2 = operand 2 guard, round and sticky bits if (OP1sex == OP2sex) { ! summation case !print "Summing^"; !font off; !print " ", (hex) OP1mhi, (hex) OP1mlo, (char) '/', (hex) tmp, " (", RNDexp, ")^"; !print "+ ", (hex) OP2mhi, (hex) OP2mlo, (char) '/', (hex) tmp2, "^"; !print " -------------^"; RNDgrs= tmp + tmp2; ! no carry possible - one of these is 0 OP1mlo = OP1mlo + OP2mlo; tmp = _UCmp(OP1mlo, OP2mlo) < 0; ! get carry tmp2 = OP1mhi + OP2mhi + tmp; tmp = (OP1mhi < 0 && OP2mhi < 0) || (OP1mhi < 0 && tmp2 >= 0) || (OP2mhi < 0 && tmp2 >= 0); OP1mhi = tmp2; ! print " ", tmp, (hex) OP1mhi, (hex) OP1mlo, (char) '/', (hex) RNDgrs, ! " (", RNDexp, ")^"; if (tmp) { @log_shift RNDgrs 1 -> tmp; RNDgrs = RNDgrs | tmp; @log_shift RNDgrs 0-1 -> RNDgrs; if (OP1mlo & 1) RNDgrs = RNDgrs | $8000; @log_shift OP1mlo 0-1 -> OP1mlo; if (OP1mhi & 1) OP1mlo = OP1mlo | $8000; @log_shift OP1mhi 0-1 -> OP1mhi; OP1mhi = OP1mhi | $8000; RNDexp = RNDexp + 1; ! print " ", (hex) OP1mhi, (hex) OP1mlo, (char) '/', (hex) RNDgrs, ! " (", RNDexp, ")^"; } !font on; } else { !print "Difference^"; !font off; !print " ", (hex) OP1mhi, (hex) OP1mlo, (char) '/', (hex) tmp, " (", RNDexp, ")^"; !print "- ", (hex) OP2mhi, (hex) OP2mlo, (char) '/', (hex) tmp2, "^"; !print " -------------^"; RNDgrs = tmp - tmp2; tmp = _UCmp(tmp, tmp2) >= 0; ! Carry tmp2 = OP1mlo - OP2mlo + tmp - 1; tmp = (OP1mlo < 0 && OP2mlo >= 0) || (OP1mlo < 0 && tmp2 >= 0) || (OP2mlo >= 0 && tmp2 >= 0); OP1mlo = tmp2; tmp2 = OP1mhi - OP2mhi + tmp - 1; tmp = (OP1mhi < 0 && OP2mhi >= 0) || (OP1mhi < 0 && tmp2 >= 0) || (OP2mhi >= 0 && tmp2 >= 0); OP1mhi = tmp2; ! print " ", 1-tmp, (hex) OP1mhi, (hex) OP1mlo, (char) '/', (hex) RNDgrs, ! " (", RNDexp, ")^"; ! If it came out negative, reverse the sign and mantissa if (~~tmp) { OP1sex = OP1sex + $8000; RNDgrs = -RNDgrs; tmp = RNDgrs == 0; tmp2 = -OP1mlo + tmp - 1; tmp = (OP1mlo >= 0 && tmp2 >= 0); OP1mlo = tmp2; OP1mhi = -OP1mhi + tmp - 1; ! print "N ", (hex) OP1mhi, (hex) OP1mlo, (char) '/', (hex) RNDgrs, ! " (", RNDexp, ")^"; } if (OP1mhi >= 0) { ! Need to normalise. Try a single bit at first, bringing the guard ! bit back into the mantissa. OP1mhi = OP1mhi + OP1mhi; if (OP1mlo < 0) OP1mhi++; OP1mlo = OP1mlo + OP1mlo; if (RNDgrs < 0) OP1mlo++; RNDgrs = RNDgrs + RNDgrs; RNDexp--; ! print " ", (hex) OP1mhi, (hex) OP1mlo, (char) '/', (hex) RNDgrs, ! " (", RNDexp, ")^"; ! If still not normalised, exponent difference must have been 0 or 1, ! so round and sticky bits are both zero. Will normalise below (in ! the same code that clears up for subnormal operands). if ((OP1mhi | OP1mlo) == 0) { ! Zero result - sign determined by rounding mode OP1sex = _rounding; if (OP1sex == 0) OP1sex = activefenv.rounding; if (OP1sex == FE_DOWNWARD) OP1sex = $8000; else OP1sex = 0; RNDexp = 0; } } } if (OP1mhi >= 0) { ! Subnormal result !if (RNDgrs ~= 0) print "[** Internal error: _fadd - RNDgrs @@126= 0 **]"; tmp = _Normalise(RNDexp, OP1mhi, OP1mlo); RNDexp = tmp-->0; OP1mhi = tmp-->1; OP1mlo = tmp-->2; ! print " ", (hex) OP1mhi, (hex) OP1mlo, (char) '/', (hex) RNDgrs, ! " (", RNDexp, ")^"; } !font on; _RoundResult(OP1sex, RNDexp, OP1mhi, OP1mlo, RNDgrs); ]; [ fmul dest OP1 OP2 rnd; _dyadic(FE_OP_MUL, _fmul, FE_FMT_S, dest, OP1, OP2, rnd); ]; [ fmulx dest OP1 OP2 rnd; _dyadic(FE_OP_MUL, _fmul, FE_FMT_X, dest, OP1, OP2, rnd); ]; [ _fmul OP1sex OP1mhi OP1mlo OP2sex OP2mhi OP2mlo tmp tmp2 r1 r2 r3 r4 t1 t2 c; ! Extract exponents tmp = OP1sex & $07FF; tmp2 = OP2sex & $07FF; ! Work out sign c = (OP1sex & $8000) + (OP2sex & $8000); if (tmp == $7FF || tmp2 == $7FF) { ! Multiplication by infinity if ((tmp2 ~= $7FF && (OP2mhi|OP2mlo) == 0) || (tmp ~= $7FF && (OP1mhi|OP1mlo) == 0)) { ! print "Infinity times zero^"; if (tmp == $7FF) tmp = InvReason_InfTimes0; else tmp = InvReason_0TimesInf; _RaiseIVO(tmp, OP1sex, OP1mhi, OP1mlo, OP2sex, OP2mhi, OP2mlo); return; } ! Return correctly signed infinity _ReturnResult(c | $07FF); return; } OP1sex = c; if (OP1mhi >= 0) { r1 = _Normalise(tmp, OP1mhi, OP1mlo); tmp = r1-->0; OP1mhi = r1-->1; OP1mlo = r1-->2; } if (OP2mhi >= 0) { r1 = _Normalise(tmp2, OP2mhi, OP2mlo); tmp2 = r1-->0; OP2mhi = r1-->1; OP2mlo = r1-->2; } if ((OP1mhi&OP2mhi) == 0) jump multzeros; OP2sex = tmp + tmp2 - 1022; !print (hex) OP1mhi, (hex) OP1mlo, " x ", (hex) OP2mhi, (hex) OP2mlo; ! new_line; ! OP1mhi * OP2mhi -> (r1,r2) _mul32(_fpscratch, OP1mhi, OP2mhi); r1 = _fpscratch-->0; r2 = _fpscratch-->1; if (OP1mlo ~= 0 && OP2mlo ~= 0) { ! OP1mlo * OP2mlo -> (r3,r4) _mul32(_fpscratch, OP1mlo, OP2mlo); r3 = _fpscratch-->0; r4 = _fpscratch-->1; } if (OP2mlo ~= 0) { ! OP1mhi * OP2mlo -> (t1, t2) _mul32(_fpscratch, OP1mhi, OP2mlo); t1 = _fpscratch-->0; t2 = _fpscratch-->1; ! Add ( 0, t1, t2, 0) ! to (r1, r2, r3, r4) r3 = r3 + t2; c = _UCmp(r3, t2) < 0; tmp = r2 + t1 + c; if ((r2 < 0 && t1 < 0) || (r2 < 0 && tmp >= 0) || (t1 < 0 && tmp >= 0)) r1++; r2 = tmp; } if (OP1mlo ~= 0) { ! OP1mlo * OP2mhi -> (t1, t2) _mul32(_fpscratch, OP1mlo, OP2mhi); t1 = _fpscratch-->0; t2 = _fpscratch-->1; ! Add ( 0, t1, t2, 0) ! to (r1, r2, r3, r4) r3 = r3 + t2; c = _UCmp(r3, t2) < 0; tmp = r2 + t1 + c; if ((r2 < 0 && t1 < 0) || (r2 < 0 && tmp >= 0) || (t1 < 0 && tmp >= 0)) r1++; r2 = tmp; } ! font off; !print " ", (hex) r1, (hex) r2, (hex) r3, (hex) r4, " (", OP2sex, ")^"; ! Make up guard, round and sticky bits: @log_shift r4 2 -> sp; @or r4 sp -> r4; @log_shift r4 0-2 -> sp; @or r3 sp -> r3; !print " ", (hex) r1, (hex) r2, "|", (hex) r3, " (", OP2sex, ")^"; if (r1 >= 0) { ! Renormalise, recovering the guard bit. r1 = r1 + r1; if (r2 < 0) ++r1; r2 = r2 + r2; if (r3 < 0) ++r2; r3 = r3 + r3; --OP2sex; ! print " ", (hex) r1, (hex) r2, "|", (hex) r3, " (", OP2sex, ")^"; } ! font on; _RoundResult(OP1sex, OP2sex, r1, r2, r3); return; .multzeros; _RoundResult(OP1sex); ]; [ fdiv dest OP1 OP2 rnd; _dyadic(FE_OP_DIV, _fdiv, FE_FMT_S, dest, OP1, OP2, rnd); ]; [ fdivx dest OP1 OP2 rnd; _dyadic(FE_OP_DIV, _fdiv, FE_FMT_X, dest, OP1, OP2, rnd); ]; [ _fdiv OP1sex OP1mhi OP1mlo OP2sex OP2mhi OP2mlo c RNDexp RNDsgn Qmhi Qmmi Qmlo bits reqbits; ! Extract exponents Qmhi = OP1sex & $07FF; Qmlo = OP2sex & $07FF; ! Work out sign RNDsgn = (OP1sex & $8000) + (OP2sex & $8000); if (Qmhi == $7FF || Qmlo == $7FF) { if (Qmlo ~= $7FF) { ! Infinity / x ! Return correctly signed infinity _ReturnResult(RNDsgn | $07FF); } else if (Qmhi ~= $7FF) { ! x / Infinity ! Return correctly signed zero _ReturnResult(RNDsgn); } else { ! Infinity / Infinity _RaiseIVO(InvReason_InfDivInf, OP1sex, OP1mhi, OP1mlo, OP2sex, OP2mhi, OP2mlo); } return; } if ((OP1mhi|OP1mlo)==0) { if ((OP2mhi|OP2mlo)==0) ! Zero / Zero _RaiseIVO(InvReason_0Div0, OP1sex, OP1mhi, OP1mlo, OP2sex, OP2mhi, OP2mlo); else ! Zero / X _ReturnResult(RNDsgn); return; } if ((OP2mhi|OP2mlo)==0) { ! X / 0 _RaiseDVZ(OP1sex, OP1mhi, OP1mlo, OP2sex); return; } if (OP1mhi >= 0) { c = _Normalise(Qmhi, OP1mhi, OP1mlo); Qmhi = c-->0; OP1mhi = c-->1; OP1mlo = c-->2; } if (OP2mhi >= 0) { c = _Normalise(Qmlo, OP2mhi, OP2mlo); Qmlo = c-->0; OP2mhi = c-->1; OP2mlo = c-->2; } ! Prospective exponent RNDexp = Qmhi - Qmlo + 1023; ! A basic long division algorithm. ! (Qmhi,Qmmi,Qmlo) will be the quotient ! (c,OP1mhi,OP1mlo) is the dividend (c using 1 bit only) if (_precision == FE_FMT_S) reqbits = 24 + 2; ! + 2 for guard + round else reqbits = 32 + 2; c=0; Qmhi=0; Qmmi=0; Qmlo=0; ! font off; for (bits = 0: bits < reqbits && (c|OP1mhi|OP1mlo): bits++) { ! print c, (hex) OP1mhi, (hex) OP1mlo, " ", (hex) OP2mhi, (hex) OP2mlo; Qmhi = Qmhi + Qmhi; if (Qmmi < 0) ++Qmhi; Qmmi = Qmmi + Qmmi; if (Qmlo < 0) ++Qmmi; Qmlo = Qmlo + Qmlo; if (c || _UCmp(OP2mhi, OP1mhi) < 0 || (OP2mhi == OP1mhi && _UCmp(OP2mlo, OP1mlo) <= 0)) { Qmlo = Qmlo | 1; c = _UCmp(OP1mlo, OP2mlo) >= 0; OP1mlo = OP1mlo - OP2mlo; OP1mhi = OP1mhi - OP2mhi + c - 1; c = 0; } ! print " ", (hex) Qmhi, (hex) Qmmi, (hex) Qmlo; new_line; c = OP1mhi < 0; OP1mhi = OP1mhi + OP1mhi; if (OP1mlo<0) ++OP1mhi; OP1mlo = OP1mlo + OP1mlo; } bits = 48 - bits; while (bits >= 16) { Qmhi = Qmmi; Qmmi = Qmlo; Qmlo = 0; bits = bits-16; } reqbits = bits-16; @log_shift Qmhi bits -> sp; @log_shift Qmmi reqbits -> sp; @or sp sp -> Qmhi; @log_shift Qmmi bits -> sp; @log_shift Qmlo reqbits -> sp; @or sp sp -> Qmmi; @log_shift Qmlo bits -> Qmlo; ! print " ", (hex) Qmhi, (hex) Qmmi, (hex) Qmlo; new_line; ! font on; if (c|OP1mhi|OP1mlo) Qmlo = Qmlo | 1; if (Qmhi>=0) { Qmhi = Qmhi + Qmhi; if (Qmmi<0) ++Qmhi; Qmmi = Qmmi + Qmmi; if (Qmlo<0) ++Qmmi; Qmlo = Qmlo + Qmlo; --RNDexp; } _RoundResult(RNDsgn, RNDexp, Qmhi, Qmmi, Qmlo); ]; [ frem dest OP1 OP2; _dyadic(FE_OP_REM, _frem, FE_FMT_S, dest, OP1, OP2); ]; [ fremx dest OP1 OP2; _dyadic(FE_OP_REM, _frem, FE_FMT_X, dest, OP1, OP2); ]; [ _frem OP1sex OP1mhi OP1mlo OP2sex OP2mhi OP2mlo OP1exp OP2exp RNDexp RNDsgn tmp tmp2 c iter REMhi; !print "_frem(", (hex)OP1sex, "|", (hex)OP1mhi, (hex)OP1mlo, ","; !print (hex)OP2sex, "|", (hex)OP2mhi, (hex)OP2mlo, ")^"; OP1exp = OP1sex & $07FF; OP2exp = OP2sex & $07FF; ! If first operand is infinity, invalid operation if (OP1exp == $7FF) { ! Inf % X _RaiseIVO(InvReason_InfRemX, OP1sex, OP1mhi, OP1mlo, OP2sex, OP2mhi, OP2mlo); return; } else if (OP1mhi >= 0) { ! If first operand is 0, return that 0. if ((OP1mhi|OP1mlo)==0 && (OP2mhi|OP2mlo)~=0) { _RoundResult(OP1sex & $8000); return; } tmp = _Normalise(OP1exp, OP1mhi, OP1mlo); OP1exp = tmp-->0; OP1mhi = tmp-->1; OP1mlo = tmp-->2; } ! If second operand is infinity, return first number unchanged if (OP2exp == $7FF) { _RoundResult(OP1sex & $8000, OP1exp, OP1mhi, OP1mlo); return; } else if (OP2mhi >= 0) { ! If second operand is 0, invalid operation if ((OP2mhi|OP2mlo)==0) { ! X % 0 _RaiseIVO(InvReason_XRem0, OP1sex, OP1mhi, OP1mlo, OP2sex, OP2mhi, OP2mlo); return; } tmp = _Normalise(OP2exp, OP2mhi, OP2mlo); OP2exp = tmp-->0; OP2mhi = tmp-->1; OP2mlo = tmp-->2; } ! Prospective exponent RNDexp = OP2exp - 1; ! Number of iterations - 1 iter = OP1exp - RNDexp; ! Isolate prospective sign, keeping a copy OP1sex = OP1sex & $8000; RNDsgn = OP1sex; if (iter < 0) { _RoundResult(RNDsgn, iter+RNDexp, OP1mhi, OP1mlo); return; } REMhi = 0; for (::) { tmp = OP2mlo - OP1mlo; c = _UCmp(OP2mlo, OP1mlo) >= 0; tmp2 = OP2mhi - OP1mhi + c - 1; if (REMhi ~= 0 || ~~((OP2mhi < 0 && OP1mhi >= 0) || (OP2mhi < 0 && tmp2 >= 0) || (OP1mhi >= 0 && tmp2 >= 0))) { OP1mlo = tmp + OP2mlo; c = _UCmp(OP1mlo, tmp) < 0; OP1mhi = tmp2 + OP2mhi + c; RNDsgn = RNDsgn + $8000; } if (--iter < 0) break; @log_shift OP1mhi 0-15 -> REMhi; OP1mhi = OP1mhi + OP1mhi; if (OP1mlo < 0) OP1mhi++; OP1mlo = OP1mlo + OP1mlo; } if ((OP1mhi|OP1mlo)==0) { ! If result is 0, should return original sign RNDsgn = OP1sex; RNDexp = 0; } else { tmp = _Normalise(RNDexp, OP1mhi, OP1mlo); RNDexp = tmp-->0; OP1mhi = tmp-->1; OP1mlo = tmp-->2; } !print "result: ", (hex) RNDexp, "|", (hex) OP1mhi, (hex) OP1mlo; new_line; _RoundResult(RNDsgn, RNDexp, OP1mhi, OP1mlo); ]; [ frnd dest OP1 rnd; _dest = dest; _precision = FE_FMT_S; _rounding = rnd; OP1 = _Promote(OP1); _frnd(OP1-->0, OP1-->1, OP1-->2); ]; [ frndx dest OP1 rnd; _dest = dest; _precision = FE_FMT_X; _rounding = rnd; _frnd(OP1-->0, OP1-->1, OP1-->2); ]; [ _frnd OP1sex OP1mhi OP1mlo RNDexp RNDsgn RNDgrs tmp c; _op = FE_OP_RND; RNDexp = OP1sex & $07FF; RNDsgn = OP1sex & $8000; if (RNDexp == $7FF) { if (OP1mhi|OP1mlo) { @test OP1mhi $4000 ?qnanorinf; _RaiseIVO(InvReason_SigNaN, OP1sex, OP1mhi, OP1mlo); return; } .qnanorinf; _ReturnResult(OP1sex, OP1mhi, OP1mlo); return; } if ((OP1mhi|OP1mlo)==0) { RNDexp = 0; jump exact; } ! tmp = position of binary point (relative to bottom of mantissa) tmp = 1023+31-RNDexp; if (tmp > 0) { ! binary point lies within or above mantissa. denormalise so that ! binary point rests at the bottom, perform normal rounding, then ! renormalise c = _Denorm(OP1mhi, OP1mlo, tmp); OP1mhi = c-->0; OP1mlo = c-->1; RNDgrs = c-->2; RNDexp = RNDexp + tmp; tmp = _precision; _precision = FE_FMT_X; c = _RoundNum(RNDsgn, RNDexp, OP1mhi, OP1mlo, RNDgrs); _precision = tmp; OP1mhi = c-->0; OP1mlo = c-->1; RNDexp = c-->2; if ((OP1mhi|OP1mlo)~=0) { c = _Normalise(RNDexp, OP1mhi, OP1mlo); RNDexp = c-->0; OP1mhi = c-->1; OP1mlo = c-->2; } else { RNDexp = 0; } } .exact; _RoundResult(RNDsgn, RNDexp, OP1mhi, OP1mlo); ]; [ fsqrt dest OP1 rnd; _dest = dest; _precision = FE_FMT_S; _rounding = rnd; OP1 = _Promote(OP1); _fsqrt(OP1-->0, OP1-->1, OP1-->2); ]; [ fsqrtx dest OP1 rnd; _dest = dest; _precision = FE_FMT_X; _rounding = rnd; _fsqrt(OP1-->0, OP1-->1, OP1-->2); ]; [ _fsqrt OP1sex OP1mhi OP1mlo RNDexp tmp tmp2 c Qhi Qlo T lastbit; _op = FE_OP_SQRT; RNDexp = OP1sex & $07FF; !print "fsqrt(", (hex) OP1sex, "|", (hex) OP1mhi, (hex) OP1mlo, ")^"; ! Normal NaN handling if (RNDexp == $7FF) { if (OP1mhi|OP1mlo) { @test OP1mhi $4000 ?qnanorinf; _RaiseIVO(InvReason_SigNaN, OP1sex, OP1mhi, OP1mlo); return; } ! sqrt(+infinity) = +infinity if (OP1sex >= 0) { .qnanorinf; _ReturnResult(OP1sex, OP1mhi, OP1mlo); return; } } else if ((OP1mhi|OP1mlo)==0) { ! sqrt(+-0) = +-0 _ReturnResult(OP1sex & $8000); return; } ! sqrt(negative) = invalid if (OP1sex < 0) { _RaiseIVO(InvReason_SqrtNeg, OP1sex, OP1mhi, OP1mlo); return; } if (OP1mhi >= 0) { tmp = _Normalise(RNDexp, OP1mhi, OP1mlo); RNDexp = tmp-->0; OP1mhi = tmp-->1; OP1mlo = tmp-->2; } ! exp+bias => (exp+2*bias)/2 = exp/2 + bias RNDexp = RNDexp + 1023; c = RNDexp & 1; @log_shift RNDexp 0-1 -> RNDexp; !font off; !print "Sqrt: m=", (hex) OP1mhi, (hex) OP1mlo; new_line; if (c == 0) { Qhi = OP1mhi - $8000; T = $1000; } else { Qhi = OP1mhi - $4000; T = $2000; } Qlo = OP1mlo; OP1mhi = $8000; OP1mlo = $0000; ! First iterations 0/1 to 13 do { !print (hex) OP1mhi, (hex) OP1mlo, " ",(hex) Qhi, (hex) Qlo, " ", (hex) T; new_line; c = Qhi < 0; Qhi = Qhi + Qhi; if (Qlo < 0) ++Qhi; Qlo = Qlo + Qlo; tmp = OP1mhi | T; if (c || _UCmp(Qhi, tmp) >= 0) { Qhi = Qhi - tmp; @log_shift T 1 -> sp; @or OP1mhi sp -> OP1mhi; } @log_shift T 0-1 -> T; } until (T==0); if ((Qhi | Qlo) == 0) jump done; if (_precision == FE_FMT_S) lastbit = $0020; else lastbit = 0; ! Iterations 14 to 23 or 14-29 T = $8000; do { !print (hex) OP1mhi, (hex) OP1mlo, " ", (hex) Qhi, (hex) Qlo, " 0000", (hex) T; new_line; c = Qhi < 0; Qhi = Qhi + Qhi; if (Qlo < 0) ++Qhi; Qlo = Qlo + Qlo; tmp = OP1mlo | T; if (c || _UCmp(Qhi, OP1mhi) > 0 || (Qhi == OP1mhi && _UCmp(Qlo, tmp) >= 0)) { Qhi = Qhi - OP1mhi; if (_UCmp(Qlo, tmp) < 0) --Qhi; Qlo = Qlo - tmp; if (T < 0) OP1mhi = OP1mhi | 1; else { @log_shift T 1 -> sp; @or OP1mlo sp -> OP1mlo; } } @log_shift T 0-1 -> T; } until (T == lastbit); !print (hex) OP1mhi, (hex) OP1mlo, " ", (hex) Qhi, (hex) Qlo; new_line; ! Iterations 30-31, if necessary if ((Qhi | Qlo) ~= 0 && lastbit == 0) { ! Manually crank out the last bit c = Qhi < 0; Qhi = Qhi + Qhi; if (Qlo < 0) ++Qhi; Qlo = Qlo + Qlo; if (c || _UCmp(Qhi, OP1mhi) > 0 || (Qhi == OP1mhi && _UCmp(Qlo, OP1mlo) > 0)) { T = 1; tmp2 = Qlo - OP1mlo - 1; c = (Qlo < 0 && OP1mlo >= 0) || (Qlo < 0 && tmp2 >= 0) || (OP1mlo >= 0 && tmp2 >= 0); Qlo = tmp2; Qhi = Qhi - OP1mhi + c - 1; OP1mlo = OP1mlo | 1; } !print (hex) OP1mhi, (hex) OP1mlo, " ", (hex) Qhi, (hex) Qlo, T; new_line; ! And the round bit c = Qhi < 0; Qhi = Qhi + Qhi; if (Qlo < 0) ++Qhi; Qlo = Qlo + Qlo + T; if (c || _UCmp(Qhi, OP1mhi) > 0 || (Qhi == OP1mhi && _UCmp(Qlo, OP1mlo) > 0)) Qhi = $8001; else Qhi = 1; Qlo = 0; !print "Round/sticky = ", (hex) Qhi; new_line; } .done; !font on; _RoundResult(OP1sex & $8000, RNDexp, OP1mhi, OP1mlo, Qhi|Qlo); ]; [ hex x y; @log_shift x 0-12 -> y; print (hexdigit) y; @log_shift x 0-8 -> y; print (hexdigit) y; @log_shift x 0-4 -> y; print (hexdigit) y; print (hexdigit) x; ]; [ hexdigit x; x = x & $F; switch (x) { 0 to 9: print (char) '0'+x; 10 to 15: print (char) 'A'-10+x; } ]; [ fraw x; print (hex) x-->0, (hex) x-->1; ]; [ frawx x; print (hex) x-->0, (char) '|', (hex) x-->1, (hex) x-->2; ]; [ fhex x; fhexx(_Promote(x)); ]; [ fhexx x exp mhi mlo tmp; exp = x-->0; mhi = x-->1; mlo = x-->2; if (exp < 0) { print (char) '-'; exp = exp - $8000; } exp = exp & $07FF; if (exp == $7FF) { if ((mhi | mlo) == 0) print "Infinity"; else { print "NaN"; if ((mhi & $4000) == 0) print (char) 'S'; mhi = mhi & $3FFF; ! Rearrange 8 extra bits of extra precision to come out at the top @log_shift mlo 8 -> sp; @log_shift sp 0-2 -> sp; @log_shift mlo 0-8 -> sp; @log_shift mhi 8 -> sp; @or sp sp -> mlo; @log_shift mhi 0-8 -> sp; @or sp sp -> mhi; print "($"; x = 0; do { @log_shift mhi 0-12 -> tmp; if (x || tmp) { x = 1; print (hexdigit) tmp; } @log_shift mhi 4 -> sp; @log_shift mlo 0-12 -> sp; @or sp sp -> mhi; @log_shift mlo 4 -> mlo; @log_shift mhi 0-12 -> tmp; } until ((mhi|mlo)==0); if (~~x) print (char) '0'; print (char) ')'; } return; } if (mhi | mlo) exp = exp - 1023 - 3; else exp = 0; @log_shift mhi 0-12 -> tmp; print (char) '$', (hexdigit) tmp; mhi = mhi & $0FFF; if (mhi | mlo) { print (char) '.'; @log_shift mhi 0-8 -> tmp; print (hexdigit) tmp; mhi = mhi & $00FF; if (mhi | mlo) { @log_shift mhi 0-4 -> tmp; print (hexdigit) tmp; mhi = mhi & $000F; if (mhi | mlo) { print (hexdigit) mhi; if (mlo) { @log_shift mlo 0-12 -> tmp; print (hexdigit) tmp; mlo = mlo & $0FFF; if (mlo) { @log_shift mlo 0-8 -> tmp; print (hexdigit) tmp; mlo = mlo & $00FF; if (mlo) { @log_shift mlo 0-4 -> tmp; print (hexdigit) tmp; mlo = mlo & $000F; if (mlo) print (hexdigit) mlo; } } } } } } print (char) 'P', exp; ]; [ fp x sxm mhi mlo tmp exp i first last dp sig wantexp; fstop(_workF0, x); sxm = _workF0-->0; mhi = _workF0-->1; mlo = _workF0-->2; !print "fp: ", (frawx) _workF0; new_line; exp = (sxm & $0FF0); @log_shift exp 0-4 -> exp; ! Turn 9 digits into ZSCII for (i=8, tmp=mlo: i>=5: i--) { _fpscratch->i = '0' + (tmp & $F); @log_shift tmp 0-4 -> tmp; } for (tmp=mhi: i>=1: i--) { _fpscratch->i = '0' + (tmp & $F); @log_shift tmp 0-4 -> tmp; } _fpscratch->0 = '0' + (sxm & $F); if (exp == $FF) { if (sxm<0) print (char) '-'; if ((mhi|mlo|(sxm & $F))==0) print "Infinity"; else { print "NaN"; if ((sxm & 7) ~= 0 || mhi ~= 0 || mlo ~= 1) { print (char) '('; if (((sxm & 7) | mhi)==0) switch (mlo) { $0000: print "SigNaN)"; return; $0004: print "MagSubInf)"; return; $0005: print "InfTimes0)"; return; $0006: print "0TimesInf)"; return; $0007: print "0Div0)"; return; $0008: print "InfDivInf)"; return; $0009: print "InfRemX)"; return; $0010: print "XRem0)"; return; $0011: print "SqrtNeg)"; return; } .asdec; _fpscratch->0 = _fpscratch->0 & $F7; for (i=0: i<8: i++) if (_fpscratch->i ~= '0') break; for (: i<9: i++) print (char) _fpscratch->i; print (char) ')'; } } return; } ! Turn exponent back into decimal @log_shift exp 0-4 -> sp; @mul sp 10 -> sp; @and exp $F -> sp; @add sp sp -> exp; if (sxm & $4000) exp = -exp; ! Find how many significant digits we have after the decimal point for (sig=8: sig>0: sig--) { if (_fpscratch->sig ~= '0') break; } ++sig; .reposition; ! Decide what the first and last digits we want are switch (activefenv.printmode) { FE_PRINT_E: first = 0; dp = 1; last = activefenv.printprecision; wantexp = true; FE_PRINT_F: if (exp >= 0) first = 0; else first = exp; dp = 1+exp; last = activefenv.printprecision + exp; wantexp = false; FE_PRINT_G: tmp = activefenv.printprecision; if (tmp<=0) tmp=1; if (exp < -4 || exp >= tmp) { first = 0; dp = 1; last = tmp-1; if (last > sig-1) last = sig-1; wantexp = true; } else { if (exp >= 0) first = 0; else first = exp; dp = 1+exp; last = tmp-1; if (last > sig-1 && last > dp-1) { if (sig > dp) last = sig-1; else last = dp-1; } wantexp = false; } } !print "first=", first, " last=", last, " sig=", sig, " dp=", dp; new_line; if (last < sig-1) { ! trailing (non-zero) digits beyond last one we're printing ! we need to round again i = _fpscratch->(last+1); sig = last+1; tmp = 0; switch (fegetround()) { FE_TONEAREST: if (i > '5' || (i == '5' && sig > (last+2)) || (i == '5' && (_fpscratch->last & 1))) tmp = 1; FE_UPWARD: tmp = 1; FE_TOWARDZERO: if (sxm < 0) tmp = 1; } if (tmp) { ! Round up - add one, looping to do carries for (i=last: i>=0: i--) { if (++(_fpscratch->i) == '9'+1) { _fpscratch->i = '0'; sig = i; } else break; } if (i<0) { ! Whoops - rounded right up _fpscratch->0 = '1'; sig = 1; exp++; } } else { ! Round down - just check trailing zeros again for (i=last: i>=1: i--) { if (_fpscratch->i == '0') sig = i; else break; } } ! Think again about what we're printing jump reposition; } ! XXX should we print -0? if (sxm < 0) print (char) '-'; for (i=first: i<=last: i++) { if (i==dp) print (char) '.'; if (i>=0 && ii; else print (char) '0'; } if (wantexp) { print (char) 'E', exp; } ]; Array _X_Ten --> $0402 $A000 $0000; ! Table look-up of powers of ten up to 10^45. [ _GetPowerOfTen dest power rnd a b c n s; n = power; s = 0; ! Halve n until it is in the range of the table while (n > 13) { @log_shift n 0-1 -> n; ++s; } ! Table of powers of ten - contains all exactly representable powers switch (n) { 0: a = $03FF; b = $8000; 1: a = $0402; b = $A000; 2: a = $0405; b = $C800; 3: a = $0408; b = $FA00; 4: a = $040C; b = $9C40; 5: a = $040F; b = $C350; 6: a = $0412; b = $F424; 7: a = $0416; b = $9896; c = $8000; 8: a = $0419; b = $BEBC; c = $2000; 9: a = $041C; b = $EE6B; c = $2800; 10: a = $0420; b = $9502; c = $F900; 11: a = $0423; b = $BA43; c = $B740; 12: a = $0426; b = $E8D4; c = $A510; 13: a = $042A; b = $9184; c = $E72A; } dest-->0 = a; dest-->1 = b; dest-->2 = c; while (s > 0) { ! Square result so far fmulx(dest, dest, dest, rnd); ! Check next bit of power --s; @sub 0 s -> sp; @log_shift power sp -> n; if (n & 1) fmulx(dest, dest, _X_Ten, rnd); } ]; [ _fstop_naninf dst sex mhi mlo rnd digits dhi dlo tmp tmp2; !print "_fstop_naninf: ", (hex)mhi, (hex)mlo; new_line; if ((mhi|mlo) ~= 0 && (mhi & $4000)==0) { ! Signalling NaN if (fpstatus & IVOE) { _workF0-->0 = sex; _workF0-->1 = mhi; _workF0-->2 = mlo; tmp = activefenv.ivo_handler; @call_vn2 tmp dst FE_FMT_P FE_OP_DEC rnd InvReason_SigNaN _workF0; } else { fpstatus = fpstatus | FE_INVALID; dst-->0 = $0FF8; dst-->1 = $0000; dst-->2 = InvReason_SigNaN; ! BCD, but okay } return; } sex = (sex & $8000) | $0FF0; if (mhi) sex = sex | $0008; mhi = mhi & $3FFF; !print "Rearranging InfNaN: ", (hex)mhi, (hex)mlo; new_line; ! Infinity/NaN ! Should trap signalling NaNs - currently caught by fstox() ! ! We do actually convert the bits of the NaN (or indeed infinity) into ! BCD. We are actually converting a single-precision NaN, and don't ! want the bottom 8 bits. So it's a conversion of 22 bits -> 7 digits. @log_shift mlo 0-8 -> sp; @log_shift mhi 8 -> sp; @or sp sp -> mlo; @log_shift mhi 0-8 -> mhi; !print "Rearranged InfNaN: ", (hex)mhi, (hex)mlo; new_line; ! Oh gawd, don't ask. This is a binary->decimal ! conversion of (mhi,mlo) -> (dhi,dlo). Each step of ! the loop divides (mhi,mlo) by ten, by using an approximation ! 1/10 = 4/5 * 1/8 ~= $0.CCCCCCCC * 1/8 ! m = $0.CCCCCCCC * i is approximated by j = i - (i>>2), k = j + (j>>4), ! l = k + (k>>8), m = l + (l>>16) ! This approvidation gives i DIV 10 <= (m>>3) <= i DIV 10 + 15, and ! we just check the remainder at the end. for (digits=0: digits<7: digits++) { tmp2 = mlo; !print "Digit ", digits; new_line; if (mhi ~= 0 || mlo < 0) { ! print "i=", (hex) mhi, (hex) mlo; new_line; @log_shift mlo 0-2 -> sp; @log_shift mhi 14 -> sp; @or sp sp -> tmp; @log_shift mhi 0-2 -> sp; @sub mhi sp -> mhi; if (_UCmp(mlo, tmp) < 0) --mhi; mlo = mlo - tmp; ! print "j=", (hex) mhi, (hex) mlo; new_line; @log_shift mlo 0-4 -> sp; @log_shift mhi 12 -> sp; @or sp sp -> tmp; mlo = mlo + tmp; @log_shift mhi 0-4 -> sp; @add mhi sp -> mhi; if (_UCmp(mlo, tmp) < 0) ++mhi; ! print "k=", (hex) mhi, (hex) mlo; new_line; @log_shift mlo 0-8 -> sp; @log_shift mhi 8 -> sp; @or sp sp -> tmp; mlo = mlo + tmp; @log_shift mhi 0-8 -> sp; @add mhi sp -> mhi; if (_UCmp(mlo, tmp) < 0) ++mhi; ! print "l=", (hex) mhi, (hex) mlo; new_line; mlo = mlo + mhi; if (_UCmp(mlo, mhi) < 0) ++mhi; ! print "m=", (hex) mhi, (hex) mlo; new_line; @log_shift mlo 0-3 -> mlo; @log_shift mhi 13 -> sp; @or mlo sp -> mlo; @log_shift mhi 0-3 -> mhi; ! print "m>>3=", (hex) mhi, (hex) mlo; new_line; @log_shift mlo 2 -> sp; @add mlo sp -> tmp; @log_shift tmp 1 -> tmp; tmp = tmp2 - tmp; ! print "remainder=", tmp; new_line; if (tmp >= 10) { tmp = tmp - 10; if (++mlo==0) ++mhi; } } else { tmp = mlo % 10; mlo = mlo / 10; } @log_shift dlo 0-4 -> dlo; @log_shift dhi 12 -> sp; @or dlo sp -> dlo; @log_shift dhi 0-4 -> dhi; @log_shift tmp 8 -> sp; @or dhi sp -> dhi; ! print "squirreled=", (hex)dhi, (hex)dlo; new_line; } !print"Converted decimal = ", (hex) dhi, (hex) dlo; new_line; dst-->0 = sex; dst-->1 = dhi; dst-->2 = dlo; ]; ! This is hard ! Binary -> decimal conversion. We only provide single-precision conversions, ! as required by the standard, using extended precision to make it work. ! ! The destination format is BCD: $SEEM $MMMM $MMMM representing ! <+/->M.MMMMMMMM x 10^(<+/->EE), M and E being BCD, top bit of S being the ! sign of the number, next bit of S being the sign of the exponent. [ fstop dst tmp rnd ! tmp = src sex mhi mlo exp arith tmp2 tmp3 tmp4 inexact digits grs c; fstox(_workF0, tmp, 1); sex = _workF0-->0; mhi = _workF0-->1; mlo = _workF0-->2; exp = sex & $07FF; if (rnd==0) rnd = activefenv.rounding; !print "fstop(", (hex) sex, "|", (hex) mhi, (hex) mlo, ")^"; if (exp == $7FF) { _fstop_naninf(dst, sex, mhi, mlo, rnd); return; } if ((mhi | mlo) == 0) { sex = sex & $8000; jump done; } ! Now have a normalised (originally single precision) number, in ! extended form. exp is in the range 1-23+(1023-127) to +254+(1023-127) ! = $36A to $47E. We now add one to the exponent (so the mantissa ! lies within [1/2 .. 1), and remove the bias. arith = exp - 1023 + 1; ! arith is now in the range [-148 .. +128]. We need to ! make it a decimal exponent. This needs a logarithm, but we'll start off ! with an approximation that can only be off by +1. ! ! We know: ! ! 2^(arith-1) <= value < 2^arith ! ! Taking base-10 logarithms: ! ! (arith-1)*log(2) <= log(value) < arith*log(2) ! ! Let log2lo and log2hi be slightly too low and high approximations to log(2). ! ! if (arith > 0): (arith-1)*log2lo <= log(value) < arith*log2hi ! if (arith <= 0): (arith-1)*log2hi <= log(value) < arith*log2lo ! ! Let D = log2hi-log2lo: ! ! if (arith > 0) ! arith*log2hi - arith*D - log2lo <= log(value) < arith*log2hi ! if (arith <= 0) ! arith*log2lo - (-arith*D) - log2hi <= log(value) < arith*log2lo ! ! Then, provided that log2lo and log2hi are such that (128*D+log2lo) <= 1 ! and (148*D+log2hi) <= 1: ! ! if (arith > 0) ! floor(arith*log2hi) - 1 <= floor(log(value)) <= floor(arith*log2hi) ! if (arith <= 0) ! floor(arith*log2lo) - 1 <= floor(log(value)) <= floor(arith*log2lo) ! ! Which gives us the desired bounds. ! ! The conditions are satisfied as long as D <= 2^(-8), but we want as ! much accuracy as we can get without overflowing a 16-bit multiplication. ! We can afford to set D to 2^(-9) giving us 8 bits of accuracy in log2, ! and 7 bits of accuracy in arith, leading to a 15-bit result. ! ! So we choose log2lo = 154 * 2^(-9) = ~ 0.30078 (log(2) = ~0.30103) ! log2hi = 155 * 2^(-9) = ~ 0.30273 if (arith > 0) tmp2 = 155; ! 2^9 * log2hi else tmp2 = 154; ! 2^9 * log2lo tmp2 = arith * tmp2; @art_shift tmp2 0-9 -> arith; !print "approximate exponent=", arith; new_line; ! Now arith-1 <= floor(log(value)) = base-10 exponent <= arith if (arith >= 0) { tmp = arith; if (arith == 0) jump expadjustdone; } else tmp = -arith; ! We now need to multiply the original value by 10^(-arith) to get ! the correct decimal mantissa. ! We'll use some FP - remember status, and disable exceptions tmp2 = fpstatus; fpstatus = 0; _GetPowerOfTen(_workF1, tmp, FE_TONEAREST); if (arith >= 0) fdivx(_workF0, _workF0, _workF1, FE_TONEAREST); else fmulx(_workF0, _workF0, _workF1, FE_TONEAREST); ! Check inexact (either in 10^tmp, or multiplication/division) inexact = fpstatus & FE_INEXACT; fpstatus = tmp2; !print "After exponent extraction: ", (frawx) _workF0; new_line; ! Get the value back sex = _workF0-->0; mhi = _workF0-->1; mlo = _workF0-->2; .expadjustdone; exp = sex & $07FF; sex = sex & $8000; digits = 9; ! Shift the mantissa so the binary point is between bits 12 and 11 of mhi ! The extra bits (not many) go into grs. exp = 1023 + 3 - exp; tmp = 16 - exp; tmp2 = -exp; @log_shift mlo tmp -> grs; @log_shift mlo tmp2 -> mlo; @log_shift mhi tmp -> sp; @or mlo sp -> mlo; @log_shift mhi tmp2 -> mhi; !print "Shifted mantissa: ", (hex) mhi, (hex) mlo, (hex) grs; new_line; ! If the mantissa is <1, decrement the arith exponent, and proceed ! to "multiply by ten", otherwise extract the first digit. exp = 0; tmp2 = 0; if (mhi & $F000) jump extract_digit; --arith; ! Stage one - three words to go, accumulating into two words do { ! First multiply by 2 mhi = mhi + mhi; if (mlo < 0) ++mhi; mlo = mlo + mlo; if (grs < 0) ++mlo; grs = grs + grs; ! Then by five - work out (mhi,mlo,grs)*4 + (mhi,mlo,grs) @log_shift grs 2 -> tmp3; @log_shift grs 0-14 -> sp; @log_shift mlo 2 -> sp; @or sp sp -> tmp4; @log_shift mlo 0-14 -> sp; @log_shift mhi 2 -> sp; @or sp sp -> tmp; grs = grs + tmp3; c = _UCmp(grs, tmp3) < 0; tmp3 = mlo + tmp4 + c; c = (mlo < 0 && tmp4 < 0) || (mlo < 0 && tmp3 >= 0) || (tmp4 < 0 && tmp3 >= 0); mlo = tmp3; mhi = mhi + tmp + c; ! print "Times 10: ", (hex) mhi, (hex) mlo, (hex) grs; new_line; .extract_digit; ! The integer part of the number is the next digit. Move it up into ! exp, and decrement the digit count. @log_shift tmp2 4 -> sp; @log_shift exp 0-12 -> sp; @or sp sp -> tmp2; @log_shift exp 4 -> sp; @log_shift mhi 0-12 -> sp; @or sp sp -> exp; mhi = mhi & $0FFF; --digits; } until (grs==0); tmp = 0; ! Second loop - two words to process in (mhi,mlo) - accumulating ! into (tmp,tmp2,exp). do { ! Multiply by 2 then 5, as before mhi = mhi + mhi; if (mlo < 0) ++mhi; mlo = mlo + mlo; @log_shift mlo 0-14 -> sp; @log_shift mlo 2 -> tmp3; mlo = mlo + tmp3; c = _UCmp(mlo, tmp3) < 0; @log_shift mhi 2 -> sp; @or sp sp -> tmp3; mhi = mhi + tmp3 + c; !print "Times 10: ", (hex) mhi, (hex) mlo; new_line; ! Extract the digit @log_shift tmp 4 -> sp; @log_shift tmp2 0-12 -> sp; @or sp sp -> tmp; @log_shift tmp2 4 -> sp; @log_shift exp 0-12 -> sp; @or sp sp -> tmp2; @log_shift exp 4 -> sp; @log_shift mhi 0-12 -> sp; @or sp sp -> exp; mhi = mhi & $0FFF; --digits; } until (digits==0); !print "Inexact=", inexact, " mhi|mlo=", (hex)mhi, (hex)mlo; new_line; ! Get final inexactitude. In practice, this doesn't work very well, ! as there are many exact cases where the two errors cancel each other ! out. inexact = inexact | mhi | mlo; @log_shift mhi 0-11 -> c; ! Round bit mhi = (mhi & $07FF) | mlo; ! Sticky bits !if (inexact || c || mhi) !{ ! if (inexact) print "Inexact "; ! if (c) print "Round "; ! if (mhi) print "Sticky "; ! new_line; !} mlo = 0; ! round up flag switch (rnd) { FE_TONEAREST: if (c) if (mhi || (exp & 1)) mlo = 1; FE_UPWARD: if (sex >= 0 && (c | mhi)) mlo = 1; FE_DOWNWARD: if (sex < 0 && (c | mhi)) mlo = 1; } if (mlo) { ! print "BCD++: ", (hex) tmp, (hex) tmp2, (hex) exp; ++exp; if ((exp & $F) == 10) { ! Need to start do a BCD carry @log_shift exp 0-1 -> c; tmp3 = c + $3333; tmp3 = (tmp3 &~ c) | (c &~ tmp3); tmp3 = tmp3 & $8888; @log_shift tmp3 0-2 -> sp; @mul sp 3 -> sp; @add exp sp -> exp; if (tmp3 & $8000) { tmp2 = tmp2 + 1; @log_shift tmp2 0-1 -> c; tmp3 = $3333 + c; tmp3 = (tmp3 &~ c) | (c &~ tmp3); tmp3 = tmp3 & $8888; @log_shift tmp3 0-2 -> sp; @mul sp 3 -> sp; @add tmp2 sp -> tmp2; if (tmp3 & $8000) { tmp = tmp + 1; @log_shift tmp 0-1 -> c; tmp3 = $3333 + c; tmp3 = (tmp3 &~ c) | (c &~ tmp3); tmp3 = tmp3 & $8888; @log_shift tmp3 0-2 -> sp; @mul sp 3 -> sp; @add tmp sp -> tmp; if (tmp & $0010) { tmp = 1; ++arith; } } } } ! print " -> ", (hex) tmp, (hex) tmp2, (hex) exp; } sex = sex | tmp; mhi = tmp2; mlo = exp; if (arith < 0) { sex = sex | $4000; arith = -arith; } tmp = 0; if (arith >= 80) { arith = arith - 80; tmp = tmp + $80; } if (arith >= 40) { arith = arith - 40; tmp = tmp + $40; } if (arith >= 20) { arith = arith - 20; tmp = tmp + $20; } if (arith >= 10) { arith = arith - 10 + $10; } tmp = tmp + arith; @log_shift tmp 4 -> tmp; sex = sex | tmp; ! print "^Final result: ", (hex) sex, (hex) mhi, (hex) mlo; new_line; .done; dst-->0 = sex; dst-->1 = mhi; dst-->2 = mlo; if (inexact) { if (fpstatus & INXE) activefenv.inx_handler(dst, FE_FMT_P, FE_OP_DEC, fegetround()); else fpstatus = fpstatus | FE_INEXACT; } ]; ! Accumulate n BCD digits from the top of src into (dest-->0,dest-->1) [ _readdigits dest src n hi lo tmp c; hi = dest-->0; lo = dest-->1; do { ! First multiply by 10 hi = hi + hi; if (lo < 0) ++hi; lo = lo + lo; @log_shift lo 0-14 -> sp; @log_shift lo 2 -> tmp; lo = lo + tmp; c = _UCmp(lo, tmp) < 0; @log_shift hi 2 -> sp; @or sp sp -> tmp; hi = hi + tmp + c; ! Then add in the new digit @log_shift src 0-12 -> tmp; lo = lo + tmp; if (_UCmp(lo, tmp) < 0) ++hi; @log_shift src 4 -> src; } until (--n == 0); dest-->0 = hi; dest-->1 = lo; ]; [ _fptos_naninf dest sex mhi mlo rnd; !print "_fptos_naninf(", (hex) sex, (hex) mhi, (hex) mlo; print ")^"; if ((mhi | mlo | (sex & $F)) == 0) { ! Infinity sex = (sex & $8000) | $7F80; jump done; } ! Check quiet / signalling NaN ! (we don't trigger signalling NaNs until converted, as we're ! supposed to supply the trap handler in extended form. Is this ! sensible?) @test sex $0008 ?~sig; @or sex $0040 -> sex; .sig; sex = (sex & $8040) | $7F80; ! Pull out the bottom 7 digits only - all we care about for NaNs _fpscratch-->0 = 0; _fpscratch-->1 = 0; @log_shift mhi 4 -> mhi; _readdigits(_fpscratch, mhi, 3); _readdigits(_fpscratch, mlo, 4); ! (mhi,mlo) = value for NaN. Could be 24 bits if unusual - knock back to ! 22, and then we're all ready mhi = _fpscratch-->0; mlo = _fpscratch-->1; !print "Read digits: ", (hex) mhi, (hex) mlo; new_line; mhi = mhi & $003F; sex = sex | mhi; if ((sex & $0040)==0) { ! We now trigger the NaN. _dest = dest; _precision = FE_FMT_S; _op = FE_OP_DEC; _rounding = rnd; _RaiseIVO(InvReason_SigNaN, sex, mhi, mlo); return; } .done; dest-->0 = sex; dest-->1 = mlo; ]; [ fptos dest src rnd sex mhi mmi mlo tmp arith; sex = src-->0; mmi = src-->1; mlo = src-->2; @log_shift sex 4 -> tmp; _fpscratch-->0 = 0; _fpscratch-->1 = 0; _readdigits(_fpscratch, tmp, 2); !print "Exponent = ", _fpscratch-->1; new_line; arith = _fpscratch-->1; if (arith > 99) { _fptos_naninf(dest, sex, mmi, mlo, rnd); return; } _fpscratch-->0 = 0; _fpscratch-->1 = sex & $F; _readdigits(_fpscratch, mmi, 4); _readdigits(_fpscratch, mlo, 4); !print "Full mantissa = ", (hex) _fpscratch-->0, (hex) _fpscratch-->1; new_line; mhi = _fpscratch-->0; mlo = _fpscratch-->1; ! Short cut for zero if ((mhi|mlo)==0) { dest-->0 = sex & $8000; dest-->1 = 0; return; } if (sex & $4000) arith = -arith; ! Adjust because we've got a 9 digit integer, not 1.8 digit decimal arith = arith - 8; ! Want to convert (mhi,mlo) into an extended precision number ! Need to normalise. We know (mhi,mlo)< 10^9 < 2^30 sex = sex & $8000; sex = sex + 1023 + 31; while (mhi >= 0) { mhi = mhi + mhi; if (mlo < 0) ++mhi; mlo = mlo + mlo; --sex; } !print (hex) sex, (char) '|', (hex) mhi, (hex) mlo; new_line; _workF0-->0 = sex; _workF0-->1 = mhi; _workF0-->2 = mlo; ! Now just need to multiply by 10^arith. |arith| <= 99, so overflow ! isn't possible (we're spared because we refuse to do binary<->decimal ! on extended precision). tmp = fpstatus; fpstatus = 0; if (arith > 0) { _GetPowerOfTen(_workF1, arith, FE_TONEAREST); fmulx(_workF0, _workF0, _workF1, FE_TONEAREST); } else if (arith < 0) { _GetPowerOfTen(_workF1, -arith, FE_TONEAREST); fdivx(_workF0, _workF0, _workF1, FE_TONEAREST); } if (fpstatus & FE_INEXACT) arith = $8000; else arith = 0; fpstatus = tmp; ! print "Extended result = ", (frawx) _workF0; new_line; ! Round back to single _precision = FE_FMT_S; _dest = dest; _rounding = rnd; _op = FE_OP_DEC; _RoundResult(_workF0-->0 & $8000, _workF0-->0 & $07FF, _workF0-->1, _workF0-->2, 0, arith); !print "Final result = ", (fraw) dest; new_line; ]; [ _strtof_inf dest str cnt len sign; --cnt; str = str + cnt; dest-->1 = $0000; if (len >= 3 && (str->1 & $DF)=='N' && (str->2 & $DF)=='F') { dest-->0 = sign | $7F80; if (len >= 8 && (str->3 & $DF)=='I' && (str->4 & $DF)=='N' && (str->5 & $DF)=='I' && (str->6 & $DF)=='T' && (str->7 & $DF)=='Y') return cnt+8; else return cnt+3; } dest-->0 = $0000; return 0; ]; [ _strtof_nan dest str cnt len sign; --cnt; str = str + cnt; if (len >= 3 && (str->1 & $DF)=='A' && (str->2 & $DF)=='N') { dest-->1 = InvReason_InitNan; if (len >= 4 && (str->3 & $DF)=='S') { dest-->0 = sign | $7F80; return cnt+4; } else { dest-->0 = sign | $7FC0; return cnt+3; } } dest-->0 = $0000; dest-->1 = $0000; return 0; ]; ! Convert ZSCII string to single. len is length of string (will ! not read beyond this). Returns number of characters consumed. [ strtof dest str len c hex sign had_dot num_ok dhi dmi dlo exp cnt exp2 negexp; ! print "strtof($", (hex) dest, ",$", (hex) str, ",", len, ")^"; if (len>=0) c = str->(cnt++); while (len>0 && c==' ') if (--len>0) c = str->(cnt++); if (len>0 && (c=='+' || c=='-')) { if (c=='-') sign = $8000; if (--len>0) c = str->(cnt++); } if (len>0) switch (c) { '$': hex = 1; if (--len>0) c = str->(cnt++); 'n','N': return _strtof_nan(dest, str, cnt, len, sign); 'i','I': return _strtof_inf(dest, str, cnt, len, sign); } while (len > 0) { !print "Got '", (char) c, "'^"; if (c=='.' && ~~had_dot) had_dot=true; else if ((c>='0' && c<='9') || (hex && ((c>='a' && c<='f') || (c>='A' && c<='F')))) { num_ok=true; if (c<='9') c=c-'0'; else if (c<='F') c=c-'A'+10; else c=c-'a'+10; if (dhi==0) { @log_shift dmi 0-12 -> dhi; @log_shift dmi 4 -> sp; @log_shift dlo 0-12 -> sp; @or sp sp -> dmi; @log_shift dlo 4 -> sp; @or sp c -> dlo; if (had_dot) --exp; } else { if (hex && c ~= '0') dlo = dlo | 1; if (~~had_dot) ++exp; } } else break; if (--len>0) c = str->(cnt++); } if (hex) exp = exp * 4; if (len>0 && num_ok && (c=='e' || c=='E' || (hex && (c=='p' || c=='P')))) { num_ok=false; if (--len>0) { c = str->(cnt++); if (c=='-' || c=='+') { if (c=='-') negexp = true; if (--len>0) c = str->(cnt++); } while (len>0 && c>='0' && c<='9') { num_ok=true; exp2 = exp2*10 + c-'0'; if (--len>0) c = str->(cnt++); } if (negexp) exp = exp-exp2; else exp = exp+exp2; } } if (len>0) cnt--; if (~~num_ok) { dest-->0=$0000; dest-->1=$0000; return 0; } !print (hex) dhi, (hex) dmi, (hex) dlo, " E", exp; new_line; if ((dhi|dmi|dlo)==0) { dest-->0=sign; ! Do we want signed zero input? Why not? dest-->1=$0000; } else { if ((dhi|dmi)==0) { dmi = dlo; dlo = 0; if (hex) exp=exp-16; else exp=exp-4; } ! print (hex) dhi, (hex) dmi, (hex) dlo, " E", exp; new_line; if (hex) { exp=exp+31+16+1023; while (dhi == 0) { dhi = dmi; dmi = dlo; dlo = 0; exp = exp - 16; } while (dhi >= 0) { dhi = dhi+dhi; if (dmi < 0) dhi++; dmi = dmi+dmi; if (dlo < 0) dmi++; dlo = dlo+dlo; exp--; } ! print (hex) dhi, (hex) dmi, (hex) dlo, " P", exp; new_line; _precision = FE_FMT_S; _dest = dest; _rounding = 0; ! FE_TONEAREST? _op = FE_OP_DEC; _RoundResult(sign, exp, dhi, dmi, dlo); } else { while (dhi==0) { @log_shift dmi 0-12 -> dhi; @log_shift dmi 4 -> sp; @log_shift dlo 0-12 -> sp; @or sp sp -> dmi; @log_shift dlo 4 -> dlo; exp--; } exp = exp+8; if (exp < -99 || exp > 99) { ! raise overflow exception (or underflow) dest-->0 = sign | $7F80; dest-->1 = 0; return cnt; } else { if (exp < 0) { exp = -exp; dhi = dhi | $4000; } dhi = dhi | sign; @div exp 10 -> sp; @log_shift sp 8 -> sp; @or dhi sp -> dhi; @mod exp 10 -> sp; @log_shift sp 4 -> sp; @or dhi sp -> dhi; _workF0-->0=dhi; _workF0-->1=dmi; _workF0-->2=dlo; ! print (frawx) _workF0; new_line; fptos(dest, _workF0); } } } return cnt; ]; #Ifndef StorageForShortName; Array StorageForShortName --> 161; #Endif; ! Initialise dest from a packed string [ finit dest str len n; @output_stream 3 StorageForShortName; print (string) str; @output_stream $FFFD; len = StorageForShortName-->0; if (len > 320) { print "[** Overlong floating-point constant **]^"; quit; } n = strtof(dest, StorageForShortName + 2, len); if (n ~= len) print "[** Floating-point constant ~", (string) str, "~ not recognized **]^"; ];