! ----------------------------------------------------------------------------
!  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 && i<sig)
      print (char) _fpscratch->i;
    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 **]^";
];