khipu/floatp.cpp

/* Definitions for the floating point types.

   This file is part of khipu.

   khipu is free software: you can redistribute it and/or modify
   it under the terms of the GNU Lesser General Public License as published by
   the Free Software Foundation; either version 3 of the License, or
   (at your option) any later version.

   This program is distributed in the hope that it will be useful,
   but WITHOUT ANY WARRANTY; without even the implied warranty of
   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
   GNU Lesser General Public License for more details.

   You should have received a copy of the GNU Lesser General Public License
   along with this program.  If not, see <https://www.gnu.org/licenses/>.  */

#include <cmath>
#include <climits>
#include <cstdio>
#include <limits>
#include <cstdlib>
#include "floatp.hpp"
#include "memory.hpp"
#include "stream.hpp"
#include "integer.hpp"
#include "bvector.hpp"
#include "io.hpp"

KP_DECLS_BEGIN

bool fnan_p (double val)
{
  return (std::fpclassify (val) == FP_NAN);
}

int finf_p (double val)
{
  return (std::fpclassify (val) != FP_INFINITE ? 0 : (val < 0 ? -1 : 1));
}

static int
fclass_aux (double val, int& sign)
{
  int ret = std::fpclassify (val);

  switch (ret)
    {
      case FP_INFINITE:
        sign = val < 0 ? -1 : 1;
        break;

      case FP_NAN:
        sign = 0;
        break;   // We don't consider -NaN

      default:
        // Assume it's a normalised number.
        ret = FP_NORMAL;
        sign = std::signbit (val);
    }

  return (ret);
}

static double
get_dbl (object obj, int& sign, int& cls)
{
  double val = as_fltobj(obj)->val;

  if (kp_unlikely (obj == FLT_PINF || obj == FLT_NINF))
    {
      cls = FP_INFINITE;
#ifdef KP_ARCH_WIDE
      sign = varobj_sign (obj);
#else
      sign = obj == FLT_NINF;
#endif
    }
  else if (kp_unlikely (obj == FLT_QNAN))
    {
      cls = FP_NAN;
      sign = 0;
    }
  else
    {
      sign = std::signbit (val);
#ifdef KP_ARCH_WIDE
      if (varobj_sign (obj))
        {
          sign ^= 1;
          val = -val;
        }
#endif
      cls = FP_NORMAL;
    }

  return (val);
}

#ifdef KP_ARCH_WIDE
#  define make_bigfloat(ptr, sign)   \
     ptrtype ((ptr), typecode::BIGFLOAT) | (sign ? SIGN_BIT : 0)

#  define F_ABS(val)   val

#else

double as_float (object obj)
{
  int sgn, cls;
  return (get_dbl (obj, sgn, cls));
}

static inline object
make_bigfloat (bigfloat *fp, int sign)
{
  if (sign)
    fp->len = -fp->len;

  return (fp->as_obj ());
}

#  define F_ABS   abs

#endif

object fltobj::make (interpreter *interp, double val)
{
  fltobj *retp = alloch<fltobj> (TYPE_SHIFT + 1);

  retp->val = val;
  retp->vo_full |= FLAGS_CONST;

  interp->alval = retp->as_obj ();
  gc_register (interp, retp);
  return (interp->alval);
}

struct fake_bigfloat : public local_varobj<bigfloat>
{
  limb_t data[LNDBL_SIZE];

  fake_bigfloat (double x)
    {
      this->local_init (this->data, KP_NELEM (this->data), 0);
      this->expo = dbltoui (fabs (x), this->data, this->len);
    }

  object get_obj (int sign)
    {
#ifdef KP_ARCH_WIDE
      return (this->as_obj () | (sign ? SIGN_BIT : 0));
#else
      if (sign)
        this->len = -this->len;
      return (this->as_obj ());
#endif
    }
};

template <typename F1, typename F2>
static result<object> ff_op (interpreter *interp, object x, object y,
                             int y_sign, F1 fn, F2 alt)
{
  int s1, s2, c1, c2, rs;
  double d1 = get_dbl (x, s1, c1);
  double d2 = get_dbl (y, s2, c2);
  double ret = KP_TRY (fn (interp, d1, d2));
  int rc = fclass_aux (ret, rs);

  if (kp_likely (rc == FP_NORMAL))
    kp_return (fltobj::make (interp, ret));
  else if ((c1 | c2) != FP_NORMAL)
    {
      if (c1 == FP_NORMAL)
        kp_return (y_sign ? neg_f (interp, y) : y);
      else if (c2 == FP_NORMAL)
        kp_return (x);

      kp_return (c1 != c2 || s1 != (s2 ^ y_sign) ? FLT_QNAN : x);
    }

  fake_bigfloat tx (d1), ty (d2);
  return (alt (interp, make_bigfloat (&tx, s1), make_bigfloat (&ty, s2)));
}

object add_ff (interpreter *interp, object x, object y)
{
  return (deref (ff_op (interp, x, y, 0,
                        [](void *, double a, double b)
                            { return (a + b); }, add_FF)));
}

object sub_ff (interpreter *interp, object x, object y)
{
  return (deref (ff_op (interp, x, y, 1,
                        [](void *, double a, double b)
                            { return (a - b); }, sub_FF)));
}

object mul_ff (interpreter *interp, object x, object y)
{
  return (deref (ff_op (interp, x, y, 0,
                        [](void *, double a, double b)
                            { return (a * b); }, mul_FF)));
}

result<object> div_ff (interpreter *interp, object x, object y)
{
  return (ff_op (interp, x, y, 0,
                 [](interpreter *ip, double a, double b) -> result<double>
                   {
                     if (b == 0)
                       return (ip->raise ("arith-error",
                                          "division by zero"));

                     return (a / b);
                   }, div_FF));
}

object neg_f (interpreter *interp, object obj)
{
  int sign, cls;
  double val = get_dbl (obj, sign, cls);

  if (cls == FP_NAN || val == 0)
    kp_return (obj);
  else if (cls == FP_INFINITE)
    kp_return (obj == FLT_PINF ? FLT_NINF : FLT_PINF);

#ifdef KP_ARCH_WIDE
  kp_return (obj ^ SIGN_BIT);
#else
  kp_return (fltobj::make (interp, -val));
#endif
}

object abs_f (interpreter *interp, object obj)
{
  if (obj == FLT_QNAN)
    kp_return (obj);

  int sign, cls;
  double val = get_dbl (obj, sign, cls);

  if (sign != 0)
#ifdef KP_ARCH_WIDE
    kp_return (obj ^ SIGN_BIT);
#else
    kp_return (fltobj::make (interp, -val));
#endif
  (void)val;
  kp_return (obj);
}

result<object> mod_ff (interpreter *interp, object x, object y)
{
  int s1, c1, s2, c2;
  double v1 = get_dbl (x, s1, c1);
  double v2 = get_dbl (y, s2, c2);

  if (v2 == 0)
    return (interp->raise ("arith-error", "modulo by zero"));
  else if (c1 == FP_NAN || c2 == FP_NAN)
    kp_return (FLT_QNAN);

  kp_return (fltobj::make (interp, fmod (v1, v2)));
}

int cmp_ff (interpreter *interp, object x, object y)
{
  int s1, c1, s2, c2;
  double v1 = get_dbl (x, s1, c1);
  double v2 = get_dbl (y, s2, c2);

  if (c1 == FP_NAN || c2 == FP_NAN)
    return (1);   // Any non-zero value will do.

  return (v1 < v2 ? -1 : (v1 > v2 ? +1 : 0));
}

bool eq_ff (interpreter *interp, object x, object y)
{
  int s1, c1, s2, c2;
  double v1 = get_dbl (x, s1, c1);
  double v2 = get_dbl (y, s2, c2);

  if ((c1 | c2) != FP_NORMAL)
    return (false);

  return (v1 == v2);
}

uint32_t hash_f (interpreter *interp, object obj)
{
  int sign, cls;
  double i, dv = get_dbl (obj, sign, cls);

  if (modf (dv, &i) != 0)
    {
      uint32_t ret = hashbuf (&dv, sizeof (dv));
      return (sign ? ~ret : ret);
    }

  int64_t iv = (int64_t)dv;
  if (fitsfixint_p (iv))
    return (hash_i (interp, fixint (iv)));

  uival uv;
  uv.qv = abs (iv);

  limb_t data[] = { uv.limbs.lo, uv.limbs.hi };
  local_varobj<bigint> bi;
  bi.local_init (data, 1 + (data[1] != 0));
#ifdef KP_ARCH_WIDE
  return (hash_I (interp, bi.as_obj () | (iv < 0 ? SIGN_BIT : 0)));
#else
  if (iv < 0)
    bi.len = -bi.len;

  return (hash_I (interp, bi.as_obj ()));
#endif
}

static result<int>
pad_writef (interpreter *interp, stream *strm, int ch,
            int width, int lv, int flags)
{
  if ((flags & (io_info::FLG_LJUST | io_info::FLG_ZERO)) || lv >= width)
    return (0);

  int rv = 0;
  char buf[256];

  lv = width - lv;
  memset (buf, ch, min (lv, (int)sizeof (buf)));
  for (; lv >= (int)sizeof (buf); lv -= sizeof (buf))
    { rv += KP_TRY (strm->write (interp, buf, sizeof (buf))); }

  rv += KP_TRY (strm->write (interp, buf, lv));
  return (rv);
}

static result<int64_t>
write_flt (interpreter *interp, stream *strm, double dbl, int width,
           int prec, int flags, int radix)
{
  uint32_t space[(DBL_MANT_DIG + 28) / 29 + 1 +
    (DBL_MAX_EXP + DBL_MANT_DIG + 28 + 8) / 9];
  uint32_t *ap, *dp, *rp, *zp;
  int e2, expo = 0, ix, jx, lx, px = 1, rv = 0;
  char buf[(DBL_MANT_DIG >> 2) + 9];
  char ebuf_st[sizeof (int) * 3];
  char *ebuf = ebuf_st + sizeof (ebuf_st), *estr;
  const char *prefix = "-0X+0X 0X-0x+0x 0x";

  if (dbl < 0)
    dbl = -dbl;
  else if (flags & io_info::FLG_SIGN)
    prefix += 3;
  else if (flags & io_info::FLG_SPACE)
    prefix += 6;
  else
    px = 0, ++prefix;

  dbl = frexpl (dbl, &e2);

  if (abs (radix) == 16)
    { // Hex-float.
      int convdig = radix < 0 ? 'A' : 'a';
      if ((dbl += dbl) != 0)
        --e2;

      if (prec != 0 && prec < (int)((DBL_DIG + 1) * .831) + 1)
        { // Round the mantissa.
          double tail = dbl;
          int q;

          for (q = prec ; ; --q)
            {
              int digit = (int)tail;
              tail -= digit;

              if (q == 0)
                {
                  if (digit % 2 != 0 ? tail >= .5 : tail > .5)
                    tail = 1 - tail;
                  else
                    tail = -tail;

                  break;
                }

              tail *= 16.l;
            }

          if (tail != 0)
            for (q = prec; q > 0; --q)
              tail *= 0.0625;

          dbl += tail;
        }

      ix = 0;   // Return value.

      if ((flags & io_info::FLG_SIGN) || px != 0)
        { ix += KP_TRY (strm->putb (interp, *prefix)); }

      int digit = (int)dbl;
      dbl -= digit;

      {
        char tmp[] = { '0', (char)(convdig - 'A' + 'X'), (char)('0' + digit) };
        ix += KP_TRY (strm->write (interp, tmp, sizeof (tmp)));
      }

      if ((flags & io_info::FLG_ALT) || dbl > 0 || prec > 0)
        {
          ix += KP_TRY (strm->putb (interp, '.'));   // XXX: Localization.

          /* The following loop must terminate, based on the above
           * assertion that FLT_RADIX is a power of 2. */
          for (; dbl > 0; --prec)
            {
              dbl *= 16.;
              digit = (int)dbl;
              dbl -= digit;

              ix += KP_TRY (strm->putb (interp, digit < 0 ?
                                        '0' : convdig - 10 + digit));
            }

          while (prec-- > 0)
            { ix += KP_TRY (strm->putb (interp, '0')); }
        }

      ix += KP_TRY (strm->putb (interp, convdig - 'A' + 'P'));
      ix += KP_TRY (strm->putb (interp, expo < 0 ? '-' : '+'));

      estr = uitostr1 (ebuf, e2, 10);
      ix += KP_TRY (strm->write (interp, estr, ebuf - estr));
      return (ix);
    }

  if ((dbl += dbl) != 0)
    {
      dbl *= 0x1p28;
      e2 -= 29;
    }

  if (prec < 0)
    prec = 6;

  {
    uint32_t *p = e2 < 0 ? space : space + KP_NELEM (space) - DBL_MANT_DIG - 1;
    ap = rp = zp = p;
  }

  do
    {
      *zp = (uint32_t)dbl;
      dbl = 1000000000 * (dbl - *zp++);
    }
  while (dbl != 0);

  while (e2 > 0)
    {
      uint32_t cy = 0;
      int sh = min (29, e2);
      for (dp = zp - 1; dp >= ap; --dp)
        {
          uint64_t twx = ((uint64_t)*dp << sh) + cy;
          *dp = (uint32_t)(twx % 1000000000);
          cy = (uint32_t)(twx / 1000000000);
        }

      if (cy != 0)
        *--ap = cy;
      for (; zp > ap && zp[-1] == 0; --zp) ;
      e2 -= sh;
    }

  while (e2 < 0)
    {
      uint32_t cy = 0;
      int sh = min (9, -e2);
      int need = 1 + (prec + DBL_MANT_DIG / 3 + 8) / 9;

      for (dp = ap; dp < zp; ++dp)
        {
          uint32_t mask = *dp & ((1u << sh) - 1);
          *dp = (*dp >> sh) + cy;
          cy = (1000000000 >> sh) * mask;
        }

      if (*ap == 0)
        ++ap;

      if (cy != 0)
        *zp++ = cy;

      uint32_t *bp = (flags & io_info::FLG_FIXED) ? rp : ap;
      if ((int)(zp - bp) > need)
        zp = bp + need;

      e2 += sh;
    }

  if (ap < zp)
    for (ix = 10, expo = 9 * (int)(rp - ap);
        *ap >= (uint32_t)ix; ix *= 10, ++expo) ;

  jx = prec - expo * !(flags & io_info::FLG_FIXED) -
    (prec != 0 && !(flags & (io_info::FLG_FIXED | io_info::FLG_EXPO)));

  if (jx < 9 * (int)(zp - rp - 1))
    {
      dp = rp + 1 + ((jx + 9 * DBL_MAX_EXP) / 9 - DBL_MAX_EXP);
      jx = (jx + 9 * DBL_MAX_EXP) % 9;

      for (ix = 10, ++jx; jx < 9; ix *= 10, ++jx) ;

      uint32_t tmp = *dp % ix;
      if (tmp != 0 || dp + 1 != zp)
        {
          double round = 2 / DBL_EPSILON;
          if (*dp % 2 != 0 || (ix == 1000000000 &&
              dp > ap && dp[-1] % 2 != 0))
            round += 2;

          double small = tmp < (uint32_t)(ix >> 1) ? 0x0.8p0 :
            tmp == (uint32_t)(ix >> 1) && dp + 1 == zp ? 0x1.0p0 : 0x1.8p0;

          if (px != 0 && *prefix == '-')
            {
              round = -round;
              small = -small;
            }

          *dp -= tmp;
          if (round + small != round)
            {
              *dp += ix;
              while (*dp > 1000000000 - 1)
                {
                  *dp-- = 0;
                  if (dp < ap)
                    *--ap = 0;

                  ++*dp;
                }

              for (ix = 10, expo = 9 * (int)(rp - ap);
                  *ap >= (uint32_t)ix; ix *= 10, ++expo) ;
            }
        }

      if (zp > dp + 1)
        zp = dp + 1;
    }

  for (; zp > ap && zp[-1] == 0; --zp) ;

  if (!(flags & (io_info::FLG_EXPO | io_info::FLG_FIXED)))
    {
      if (prec == 0)
        prec = 1;

      if (prec > expo && expo >= -4)
        {
          flags |= io_info::FLG_FIXED;
          prec -= expo + 1;
        }
      else
        {
          flags |= io_info::FLG_EXPO;
          --prec;
        }

      if (!(flags & io_info::FLG_ALT))
        {
          if (zp > ap && zp[-1] != 0)
            for (ix = 10, jx = 0; zp[-1] % ix == 0; ix *= 10, ++jx) ;
          else
            jx = 9;

          int uoff = 9 * ((int)(zp - rp) - 1) - jx;
          if (flags & io_info::FLG_EXPO)
            uoff += expo;

          prec = min (prec, max (0, uoff));
        }
    }

  lx = prec || (flags & io_info::FLG_ALT);
  if (prec > INT_MAX - 1 - lx)
    return (-1);

  lx += prec + 1;

  if (flags & io_info::FLG_FIXED)
    {
      if (expo > INT_MAX - lx)
        return (-1);
      else if (expo > 0)
        lx += expo;
    }
  else
    {
      estr = uitostr1 (ebuf, abs (expo), 10);
      while (ebuf - estr < 2)
        *--estr = '0';

      *--estr = expo < 0 ? '-' : '+';
      *--estr = radix < 0 ? 'E' : 'e';
      lx += ebuf - estr;
    }

  if (lx > INT_MAX - px)
    return (-1);

  rv += KP_TRY (pad_writef (interp, strm, ' ', width, px + lx, flags));
  rv += KP_TRY (strm->write (interp, prefix, px));
  rv += KP_TRY (pad_writef (interp, strm, '0', width,
                            px + lx, flags ^ io_info::FLG_ZERO));

  if (flags & io_info::FLG_FIXED)
    {
      if (ap > rp)
        ap = rp;

      for (dp = ap; dp <= rp; ++dp)
        {
          char *s = uitostr1 (buf + 9, *dp, 10);
          if (dp != ap)
            while (s > buf)
              *--s = '0';
          else if (s == buf + 9)
            *--s = '0';

          rv += KP_TRY (strm->write (interp, s, buf + 9 - s));
        }

      if (prec || (flags & io_info::FLG_ALT))
        { rv += KP_TRY (strm->putb (interp, '.')); }   // XXX: Localization.

      for (; dp < zp && prec > 0; ++dp, prec -= 9)
        {
          char *s = uitostr1 (buf + 9, *dp, 10);
          while (s > buf)
            *--s = '0';

          rv += KP_TRY (strm->write (interp, s, min (9, prec)));
        }

      rv += KP_TRY (pad_writef (interp, strm, '0', prec + 9, 9, 0));
    }
  else
    {
      if (zp <= ap)
        zp = ap + 1;

      for (dp = ap; dp < zp && prec >= 0; ++dp)
        {
          char *s = uitostr1 (buf + 9, *dp, 10);
          if (s == buf + 9)
            *--s = '0';
          if (dp != ap)
            while (s > buf)
              *--s = '0';
          else
            {
              rv += KP_TRY (strm->putb (interp, *s++));
              if (prec > 0 || (flags & io_info::FLG_ALT))
                { rv += KP_TRY (strm->putb (interp, '.')); }
            }

          rv += KP_TRY (strm->write (interp, s,
                                     min ((int)(buf + 9 - s), prec)));
          prec -= (int)(buf + 9 - s);
        }

      rv += KP_TRY (pad_writef (interp, strm, '0', prec + 18, 18, 0));
      rv += KP_TRY (strm->write (interp, estr, ebuf - estr));
    }

  rv += KP_TRY (pad_writef (interp, strm, ' ', width,
                            px + lx, flags ^ io_info::FLG_LJUST));
  return (rv);
}

result<int64_t> write_f (interpreter *interp, stream *strm,
                         object obj, io_info& info)
{
  int sign, cls;
  double v = get_dbl (obj, sign, cls);

  if (cls == FP_INFINITE)
    {
      char buf[] = { sign ? '-' : '+', 'I', 'N', 'F' };
      return (strm->write (interp, buf, 4));
    }
  else if (cls != FP_NORMAL)
    return (strm->write (interp, "NaN", 3));

  return (write_flt (interp, strm, v, info.width, info.prec,
                     info.flags, info.radix));
}

result<int64_t> pack_f (interpreter *interp, stream *strm,
                        object obj, pack_info&)
{
  int sign, cls;
  double v = get_dbl (obj, sign, cls);
  return (strm->write (interp, &v));
}

result<object> unpack_f (interpreter *interp, stream *strm,
                         pack_info& info, bool)
{
  double v;
  bool rv = KP_TRY (strm->sread (interp, &v));
  if (rv)
    kp_return (fltobj::make (interp, v));
  return (info.error ("failed to read floating point value"));
}

// Big float implementation.

static_assert ((alignof (bigfloat) % alignof (limb_t)) == 0,
               "invalid alignment for big floats");

bigfloat* bigfloat::alloc_raw (int len)
{
  bigfloat *retp = (bigfloat *)alloch (sizeof (*retp) + sizeof (limb_t) * len,
                                       typecode::BIGFLOAT, TYPE_SHIFT + 1);

  retp->len = len;
  retp->expo = 0;
  retp->data = (limb_t *)&retp[1];
  retp->vo_full = FLAGS_CONST;

  return (retp);
}

object alloc_bigfloat (interpreter *interp, int len)
{
  auto retp = bigfloat::alloc_raw (len);
  interp->alval = retp->as_obj ();
  gc_register (interp, retp, sizeof (*retp) + len * sizeof (limb_t));
  return (interp->alval);
}

static bigfloat*
get_bigfloat (object obj, int& sign)
{
#ifdef KP_ARCH_WIDE
  sign = (obj & SIGN_BIT) / SIGN_BIT;
#else
  sign = as_bigfloat(obj)->len < 0;
#endif
  return (as_bigfloat (obj));
}

static_assert (sizeof (bigfloat) >= sizeof (fltobj),
               "big floats must be larger or equal in size to floats");

static object
ret_F (interpreter *interp, bigfloat *src, int sign)
{
  if (src->len == 0)
    {
      xfree (src);
      return (FLT_ZERO);
    }

  int b2exp = (src->expo - src->len) * LIMB_BITS;
  if (b2exp >= DBL_MIN_EXP && b2exp <= DBL_MAX_EXP)
    { // May fit in a regular double.
      double dbl = uitodbl (src->data, src->len, src->expo);
      if (finf_p (dbl))
        goto do_bigfloat;

      /* Fits in a regular double - Transmute the type.
       * This wastes some memory, but is preferable to
       * allocating a new object. */

      fltobj tmp;
      tmp.gc_link = src->gc_link;
      tmp.vo_full = src->vo_full;
      tmp.vo_type = typecode::FLOAT;
      tmp.val = sign ? -dbl : dbl;

      memcpy (src, &tmp, sizeof (tmp));
      interp->retval = src->as_obj ();
    }
  else
  do_bigfloat:
    interp->retval = make_bigfloat (src, sign);

  gc_register (interp, src, sizeof (*src) + src->len * sizeof (limb_t));
  return (interp->retval);
}

static bigfloat*
bfloat_add (interpreter *interp, const bigfloat *x, const bigfloat *y)
{ 
  if (x->expo < y->expo)
    swap (x, y);

  int ediff = x->expo - y->expo, xl = F_ABS (x->len), yl = F_ABS (y->len);
  int prec = max (xl, yl);
  if (ediff >= prec)
    // Y completely cancelled.
    return ((bigfloat *)x);

  auto ret = bigfloat::alloc_raw (prec + 1);
  ret->expo = x->expo;

  int size;
  limb_t cy;

  if (!(xl > ediff))
    {
      size = yl + ediff - xl;
      memcpy (ret->data, y->data, yl * sizeof (limb_t));
      memset (ret->data + yl, 0, (ediff - xl) * sizeof (limb_t));
      memcpy (ret->data + size, x->data, xl * sizeof (limb_t));
      cy = 0, ret->len = size + xl;
    }
  else if (yl + ediff <= xl)
    {
      size = xl - ediff - yl;
      memcpy (ret->data, x->data, size * sizeof (limb_t));
      cy = uiadd (ret->data + size, x->data + size, xl - size, y->data, yl);
      ret->len = yl;
    }
  else
    {
      size = yl + ediff - xl;
      memcpy (ret->data, y->data, size * sizeof (limb_t));
      cy = uiadd (ret->data + size, x->data, xl, y->data + size, xl - ediff);
      ret->len = yl + ediff;
    }

  ret->data[ret->len] = cy;
  ret->len += cy, ret->expo += cy;
  return (ret);
}

static bigfloat*
bfloat_sub (interpreter *interp, const bigfloat *x,
            const bigfloat *y, int *negp)
{
  int neg = 0;

  if (x->expo < y->expo)
    {
      swap (x, y);
      neg = 1;
    }

  int expo = x->expo, ediff = x->expo - y->expo;
  int xl = F_ABS (x->len), yl = F_ABS (y->len), len = max (xl, yl);
  limb_t *xp = x->data, *yp = y->data;
  auto ret = bigfloat::alloc_raw (len);
  limb_t cy;

  if (ediff <= 1)
    { /* The numbers are very close. Scan them and
       * ignore the limbs that compare equal. */
      if (ediff == 0)
        {
          if (xp[xl - 1] == yp[yl - 1])
            do
              {
                --yl, --expo;
                if (--xl == 0)
                  {
                    neg ^= 1;
                  cancel:
                    for (; yl && !yp[yl - 1]; --yl, --expo) ;
                    memcpy (ret->data, y->data, yl * sizeof (limb_t));
                    ret->len = yl;
                    goto done;
                  }

                if (yl == 0)
                  {
                    yp = xp;
                    yl = xl;
                    goto cancel;
                  }
              }
            while (xp[xl - 1] == yp[yl - 1]);

          if (xp[xl - 1] < yp[yl - 1])
            {
              swap (xp, yp);
              swap (xl, yl);
              neg ^= 1;
            }

          if (xp[xl - 1] != yp[yl - 1] + 1)
            goto main;
          
          --xl, --yl, --expo;
        }
      else   // (ediff == 1)
        {
          if (xp[xl - 1] != 1 || yp[yl - 1] != ~(limb_t)0 ||
              (xl >= 2 && xp[xl - 2] != 0))
            goto main;

          --xl, --expo;
        }

      for (; yl && xl && !xp[xl - 1] &&
          yp[yl - 1] == ~(limb_t)0; --xl, --yl, --expo) ;

      if (xl == 0)
        for (; yl && yp[yl - 1] == ~(limb_t)0; --yl, --expo) ;
      
      if (yl == 0)
        {
          memcpy (ret->data, xp, (ret->len = xl) * sizeof (limb_t));
          cy = 0;
        }
      else if (xl == 0)
        {
          for (int i = 0; i < yl; ++i)
            ret->data[i] = ~yp[i];
          cy = 1 - uiadd1 (ret->data, ret->data, yl, 1);
          ret->len = yl;
        }
      else if (xl >= yl)
        {
          int size = xl - yl;
          memcpy (ret->data, xp, size * sizeof (limb_t));
          cy = uisubn (ret->data + size, xp + size, yp, yl);
          ret->len = xl;
        }
      else
        {
          int size = yl - xl;
          for (int i = 0; i < size; ++i)
            ret->data[i] = ~yp[i];

          cy  = uisubn (ret->data + size, xp, yp + size, xl);
          cy += uisub1 (ret->data + size, ret->data + size, xl, 1);
          cy -= uiadd1 (ret->data, ret->data, yl, 1);
          ret->len = yl;
        }

      if (cy == 0)
        {
          ret->data[ret->len++] = 1;
          ++expo;
          goto normalize;
        }
    }
  
main:
  if (ediff >= len)
    {
      memcpy (ret->data, xp, xl);
      ret->len = xl;
      expo = ((const bigfloat *)
        ((char *)xp - sizeof (bigfloat)))->expo;
      goto done;
    }

  // Search for the least significant non-zero limbs in X and Y.
  
  while (true)
    {
      if (yl == 0)
        {
          memcpy (ret->data, xp, xl * sizeof (limb_t));
          ret->len = xl;
          goto done;
        }
      
      if (*yp != 0)
        break;
      
      ++yp, --yl;
    }

  while (true)
    {
      if (xl == 0)
        {
          memcpy (ret->data, yp, yl * sizeof (limb_t));
          ret->len = yl;
          neg ^= 1;
          goto done;
        }
      
      if (*xp != 0)
        break;
      
      ++xp, --xl;
    }

  if (xl > ediff)
    { // Partial overlapping.
      if (!ediff)
        { /* Compare leading limbs to determine whether to
           * do X - Y or Y - X. */
          if (xl >= yl)
            {
              int size = xl - yl;
              memcpy (ret->data, xp, size * sizeof (limb_t));
              uisubn (ret->data + size, xp + size, yp, yl);
              ret->len = xl;
            }
          else
            {
              int size = yl - xl;
              *ret->data = -*yp;
              for (int i = 1; i < size; ++i)
                ret->data[i] = ~yp[i];
              
              uisubn (ret->data + size, xp, yp + size, xl);
              uisub1 (ret->data + size, ret->data + size, xl, 1);
              ret->len = yl;
            }
        }
      else if (yl + ediff <= xl)
        {
          int size = xl - ediff - yl;
          memcpy (ret->data, xp, size * sizeof (limb_t));
          uisub (ret->data + size, xp + size, xl - size, yp, yl);
          ret->len = xl;
        }
      else
        {
          int size = yl + ediff - xl;
          *ret->data = -*yp;
          for (int i = 1; i < size; ++i)
            ret->data[i] = ~yp[i];

          uisub (ret->data + size, xp, xl, yp + size, xl - ediff);
          uisub1 (ret->data + size, ret->data + size, xl, 1);
          ret->len = yl + ediff;
        }
    }
  else
    { // No overlapping - Trivial copy.
      int size = yl + ediff - xl;
      *ret->data = -*yp;
      for (int i = 1; i < yl; ++i)
        ret->data[i] = ~yp[i];
      memset (&ret->data[yl], 0xff,
              (size - yl) * sizeof (limb_t));

      uisub1 (ret->data + size, xp, xl, 1);
      ret->len = size + xl;
    }

normalize:
  for (; ret->len && !ret->data[ret->len - 1]; --ret->len, --expo) ;

done:
  ret->expo = !ret->len ? 0 : expo;
  if (negp != nullptr) *negp = neg;
  return (ret);
}

static bigfloat*
bfloat_mul (interpreter *interp, const bigfloat *x, const bigfloat *y)
{
  int xl = F_ABS (x->len), yl = F_ABS (y->len), tsize = xl + yl;
  auto ret = bigfloat::alloc_raw (tsize);

  limb_t cy = x == y ? uisqr (interp, ret->data, x->data, xl) :
                       uimul (interp, ret->data, x->data, xl, y->data, yl);

  cy = !cy;
  ret->len -= cy;
  ret->expo = x->expo + y->expo - cy;
  return (ret);
}

static bigfloat*
bfloat_div (interpreter *interp, const bigfloat *x, const bigfloat *y)
{
  int xl = F_ABS (x->len), yl = F_ABS (y->len), rsize = xl - yl;
  limb_t *tmp, *xp = x->data;
  tmp_allocator ta { interp };

  if (rsize <= 0)
    { // Copy and extend X.
      uint32_t padding = rsize == 0 ? 1 : (rsize = -rsize) * 2;
      xp = (limb_t *)ta.alloc ((xl * 2 + padding) * sizeof (*xp));
      memset (xp, 0, padding * sizeof (limb_t));
      memcpy (xp + padding, x->data, xl * sizeof (limb_t));
      xl += padding, tmp = xp + xl;
    }
  else
    tmp = (limb_t *)ta.alloc ((xl + 1) * sizeof (*tmp));

  rsize = uidiv (interp, xp, xl, y->data, yl, tmp, false);

  auto ret = bigfloat::alloc_raw (rsize);
  memcpy (ret->data, tmp + yl, ret->len * sizeof (limb_t));
  ret->expo = x->expo - y->expo + 1 -
              (xl - yl + 1 != ret->len);   // i.e: msd equals zero.

  return (ret);
}

object add_FF (interpreter *interp, object x, object y)
{
  int s1, s2;
  auto v1 = get_bigfloat (x, s1), v2 = get_bigfloat (y, s2);
  bigfloat *ret;

  if (!(s1 ^ s2))
    {
      ret = bfloat_add (interp, v1, v2);
      if (ret == v1)
        kp_return (x);
      else if (ret == v2)
        kp_return (y);
    }
  else
    ret = s1 ? bfloat_sub (interp, v2, v1, &s1) :
               bfloat_sub (interp, v1, v2, &s1);

  return (ret_F (interp, ret, s1));
}

object add_fF (interpreter *interp, object x, object y)
{
  int s1, c1, s2;
  double v1 = get_dbl (x, s1, c1);
  const auto v2 = get_bigfloat (y, s2);

  if (c1 != FP_NORMAL)
    kp_return (x);   // NaN or +/- Inf.

  fake_bigfloat tmp (v1);
  bigfloat *ret;

  if (!(s1 ^ s2))
    {
      ret = bfloat_add (interp, &tmp, v2);
      if (ret == (bigfloat *)&tmp)
        kp_return (x);
      else if (ret == v2)
        kp_return (y);
    }
  else
    ret = s1 ? bfloat_sub (interp, v2, &tmp, &s1) :
      bfloat_sub (interp, &tmp, v2, &s1);

  return (ret_F (interp, ret, s1));
}

object sub_FF (interpreter *interp, object x, object y)
{
#ifdef KP_ARCH_WIDE
  return (add_FF (interp, x, y ^ SIGN_BIT));
#else
  local_varobj<bigfloat> tmp;
  bigfloat *yp = as_bigfloat (y);

  tmp.local_init (yp->data, -yp->len, yp->expo);
  return (add_FF (interp, x, tmp.as_obj ()));
#endif
}

object sub_fF (interpreter *interp, object x, object y)
{
#ifdef KP_ARCH_WIDE
  return (add_fF (interp, x, y ^ SIGN_BIT));
#else
  local_varobj<bigfloat> tmp;
  bigfloat *yp = as_bigfloat (y);

  tmp.local_init (yp->data, -yp->len, yp->expo);
  return (add_fF (interp, x, tmp.as_obj ()));
#endif
}

object neg_F (interpreter *interp, object obj)
{
#ifdef KP_ARCH_WIDE
  kp_return (obj ^ SIGN_BIT);
#else
  const bigfloat *srcp = as_bigfloat (obj);
  int sl = F_ABS (srcp->len);
  bigfloat *retp = as_bigfloat (alloc_bigfloat (interp, sl));

  memcpy (retp->data, srcp->data, sl * sizeof (limb_t));
  retp->expo = srcp->expo;
  retp->len = -sl;
  kp_return (retp->as_obj ());
#endif
}

object mul_FF (interpreter *interp, object x, object y)
{
  int s1, s2;
  const auto v1 = get_bigfloat (x, s1), v2 = get_bigfloat (y, s2);

  return (ret_F (interp, bfloat_mul (interp, v1, v2), s1 ^ s2));
}

object mul_fF (interpreter *interp, object x, object y)
{
  int s1, c1, s2;
  double v1 = get_dbl (x, s1, c1);

  if (c1 != FP_NORMAL)
    kp_return (x);   // NaN or +/- Inf.

  fake_bigfloat tmp (v1);
  const auto v2 = get_bigfloat (y, s2);
  return (ret_F (interp, bfloat_mul (interp, &tmp, v2), s1 ^ s2));
}

object div_FF (interpreter *interp, object x, object y)
{
  int s1, s2;
  const auto v1 = get_bigfloat (x, s1), v2 = get_bigfloat (y, s2);
  return (ret_F (interp, bfloat_div (interp, v1, v2), s1 ^ s2));
}

#define DIV_BODY(left, right, op, alt, zero_x)   \
  int s1, s2, cl;   \
  double v1 = get_dbl (left, s1, cl);   \
  const auto v2 = get_bigfloat (right, s2);   \
     \
  if (cl != FP_NORMAL)   \
    kp_return (alt);   \
  else if (v1 == 0)   \
	kp_return (zero_x);   \
     \
  fake_bigfloat tmp (v1);   \
  return (ret_F (interp, op, s1 ^ s2))

object div_fF (interpreter *interp, object x, object y)
{
  DIV_BODY (x, y, bfloat_div (interp, &tmp, v2), x, FLT_ZERO);
}

object div_Ff (interpreter *interp, object x, object y)
{
  DIV_BODY (y, x, bfloat_div (interp, v2, &tmp),
            cl == FP_INFINITE ? FLT_ZERO : x, FLT_PINF);
}

#undef DIV_BODY

static bigfloat*
bfloat_mod (interpreter *interp, const bigfloat *x, const bigfloat *y)
{
  bigfloat *p, *q, *r;

  p = bfloat_div (interp, x, y);   // p := x / y

  // Truncate P in place.
  if (p->expo <= 0)
    p->expo = p->len = 0;
  else
    {
      int len = min (p->len, p->expo);
      limb_t *lp = p->data + p->len - len;

      if (lp != p->data)
        memmove (p->data, lp, (p->len - len) * sizeof (*lp));
      
      p->len = len;
    }

  q = bfloat_mul (interp, p, y);   // q := trunc (p) * y
  r = bfloat_sub (interp, x, q, nullptr);   // r := x - q
  xfree (p);
  xfree (q);
  return (r);
}

object mod_FF (interpreter *interp, object x, object y)
{
  int s1, s2;
  const auto v1 = get_bigfloat (x, s1), v2 = get_bigfloat (y, s2);
  return (ret_F (interp, bfloat_mod (interp, v1, v2), s1));
}

object mod_fF (interpreter *interp, object x, object y)
{
  int s1, c1, s2;
  double v1 = get_dbl (x, s1, c1);

  if (c1 != FP_NORMAL)
    kp_return (FLT_QNAN);

  fake_bigfloat tmp (v1);
  const auto v2 = get_bigfloat (y, s2);
  return (ret_F (interp, bfloat_mod (interp, &tmp, v2), s1));
}

object mod_Ff (interpreter *interp, object x, object y)
{
  int s1, c2, s2;
  double v2 = get_dbl (y, s2, c2);

  if (c2 != FP_NORMAL)
    kp_return (FLT_QNAN);

  const auto v1 = get_bigfloat (x, s1);
  fake_bigfloat tmp (v2);
  return (ret_F (interp, bfloat_mod (interp, v1, &tmp), s1));
}

static int
bfloat_cmp (const bigfloat *x, int s1, const bigfloat *y, int s2)
{
  int ret = 0;

  if (s1 != s2)
    ret = 1;
  else if (!(ret = x->expo - y->expo))
    {
      limb_t *xp = x->data, *yp = y->data;
      int xl, yl;

      for (; *xp == 0; ++xp) ;
      for (; *yp == 0; ++yp) ;
      xl = (int)((x->data + x->len) - xp);
      yl = (int)((y->data + y->len) - yp);
      
      if (xl > yl)
        ret = uicmpn (xp + xl - yl, yp, yl);
      else if (yl > xl)
        ret = -uicmpn (xp, yp + yl - xl, xl);
      else
        ret = uicmpn (xp, yp, yl);
    }

  return (s1 ? -ret : ret);
}

int cmp_FF (interpreter *interp, object x, object y)
{
  int s1, s2;
  const auto v1 = get_bigfloat (x, s1), v2 = get_bigfloat (y, s2);
  return (bfloat_cmp (v1, s1, v2, s2));
}

int cmp_fF (interpreter *interp, object x, object y)
{
  int s1, c1, s2;
  double v1 = get_dbl (x, s1, c1);

  if (c1 != FP_NORMAL)
    return (s1);   // i.e: +/-1 for Inf, +1 for NaN.

  const bigfloat *v2 = get_bigfloat (y, s2);
  fake_bigfloat tmp (v1);
  return (bfloat_cmp (&tmp, s1, v2, s2));
}

bool eq_FF (interpreter *interp, object x, object y)
{
  int s1, s2;
  const auto v1 = get_bigfloat (x, s1), v2 = get_bigfloat (y, s2);

  // Exponents and signs must match.
  if (s1 != s2 || v1->expo != v2->expo)
    return (false);
  else if (v1->data == v2->data)
    return (true);

  int minlen, maxlen, i, len, diff;
  int xl = F_ABS (v1->len), yl = F_ABS (v2->len);
  const limb_t *xp = v1->data + xl, *yp = v2->data + yl;
  // Test for most significant limb position.
  if ((*(yp - 1) >> (LIMB_BITS - uiclz (*(xp - 1) - 1))) != 1)
    return (0);

  if (xl < yl)
    minlen = xl, maxlen = yl;
  else
    minlen = yl, maxlen = xl;

  // Compare the most significant limbs.
  xp += minlen, yp += minlen;
  for (i = minlen - 1; i > 0; --i)
    if (xp[i] != yp[i])
      return (0);

  if (0 != (len = maxlen - minlen))
    { /* Either X or Y has been completely examined.
       * Test that the other is made of nil limbs. */
      if (*xp != *yp)
        return (0);
      
      const limb_t *p = (xl > yl ? xp : yp) - len;
      for (i = len -1 ; i > 0; --i)
        if (p[i] != 0)
          return (0);
      
      diff = *p;
    }
  else
    diff = *xp ^ *yp;

  return (diff == 0);
}

static object
bfloat_ceilfloor (interpreter *interp, object obj, int dir)
{
  int osign;
  const auto x = get_bigfloat (obj, osign);
  int sign = osign ? -1 : 1;

  if (x->expo <= 0)
    kp_return ((sign ^ dir) < 0 ? FLT_ZERO : fltobj::make (interp, dir));

  int len = min (F_ABS (x->len), x->expo);
  auto ret = bigfloat::alloc_raw (len);
  limb_t *xp = x->data + x->len - len;
      
  ret->expo = x->expo;
  if ((sign ^ dir) >= 0)
    {
      const limb_t *tp;
      for (tp = x->data; tp != xp; ++tp)
        if (*tp != 0)
          {
            if (uiadd1 (ret->data, xp, len, 1))
              {
                ret->len = *ret->data = 1;
                ++ret->expo;
              }
                
            return (ret_F (interp, ret, osign));
          }
    }
      
  memcpy (ret->data, xp, len * sizeof (limb_t));
  return (ret_F (interp, ret, osign));
}

object floor_F (interpreter *interp, object obj)
{
  return (bfloat_ceilfloor (interp, obj, -1));
}

object ceil_F (interpreter *interp, object obj)
{
  return (bfloat_ceilfloor (interp, obj, +1));
}

object trunc_F (interpreter *interp, object obj)
{
  int sign;
  const auto x = get_bigfloat (obj, sign);

  if (x->expo <= 0)
    kp_return (FLT_ZERO);

  int xl = F_ABS (x->len), len = min (xl, x->expo); 
  auto ret = bigfloat::alloc_raw (len);
  memcpy (ret->data, x->data + xl - len, len * sizeof (limb_t));
  ret->expo = x->expo;
  return (ret_F (interp, ret, sign));
}

uint32_t hash_F (interpreter *interp, object obj)
{
  int sign;
  const auto lp = get_bigfloat (obj, sign);
  int fl = F_ABS (lp->len);

  if (fl == lp->expo)
    {
      local_varobj<bigint> bi;
      bi.local_init (lp->data, lp->len);

#ifdef KP_ARCH_WIDE
      return (hash_I (interp, bi.as_obj () | (obj & SIGN_BIT)));
#else
      return (hash_I (interp, bi.as_obj ()));
#endif
    }

  uint32_t ret = mix_hash (2 * fl, lp->expo);

  for (int i = 0; i < fl; ++i)
    ret = mix_hash (ret, lp->data[i]);

  return (sign ? ~ret :ret);
}

result<int64_t> write_F (interpreter *interp, stream *strm,
                         object obj, io_info& info)
{
  const auto lp = as_bigfloat (obj);
  int expo = lp->expo, xl = F_ABS (lp->len), ret = 0;
  int size = uibsize (info.radix, lp->data[xl - 1], xl);
  tmp_allocator ta { interp };
  char *ptr = (char *)ta.alloc (size);

  size = lftostr (interp, ptr, &expo, lp->data, xl, size, info.radix);

#ifdef KP_ARCH_WIDE
  if (obj & SIGN_BIT)
#else
  if (lp->len < 0)
#endif
    { ret += KP_TRY (strm->putb (interp, '-')); }

  // Write the decimal separator, followed by the mantissa.
  ret += KP_TRY (strm->write (interp, "0.", 2));
  ret += KP_TRY (strm->write (interp, ptr, size));

  // Write the exponent.
  size = sprintf (ptr, "%c%d", info.radix <= 10 ? 'e' : '@', expo);
  ret += KP_TRY (strm->write (interp, ptr, size));
  return (ret);
}

result<int64_t> pack_F (interpreter *interp, stream *strm,
                        object obj, pack_info&)
{
  const bigfloat *lp = as_bigfloat (obj);
  int expo = lp->expo, xl = lp->len;
  int64_t ret = 0;

#ifdef KP_ARCH_WIDE
  if (obj & SIGN_BIT)
    xl = -xl;
#endif

  ret += KP_TRY (strm->write (interp, &expo));
  ret += KP_TRY (strm->write (interp, &xl));
  ret += KP_TRY (strm->write (interp, lp->data,
                              F_ABS (lp->len) * sizeof (*lp->data)));
  return (ret);
}

result<object> unpack_F (interpreter *interp, stream *strm,
                         pack_info& info, bool save)
{
  int vals[2];
  {
    auto tmp = KP_TRY (strm->read (interp, vals, sizeof (vals)));
    if ((size_t)tmp != sizeof (vals))
      return (info.error ("failed to read floating point value"));
  }

  int len = abs (vals[1]);
  bigfloat *lp = as_bigfloat (alloc_bigfloat (interp, len));
  auto tmp = KP_TRY (strm->read (interp, lp->data, len * sizeof (*lp->data)));
  if ((size_t)tmp != (size_t)(len * sizeof (*lp->data)))
    return (info.error ("invalid number of limbs read"));
  else if (save)
    KP_VTRY (info.add_mapping (interp, *info.offset, lp->as_obj ()));

  lp->expo = vals[0];
  kp_return (make_bigfloat (lp, vals[1] < 0));
}

// Mixed operations.

static inline object
int2flt (object i, local_varobj<fltobj>& out)
{
  out.vo_full = 0;
  out.vo_type = typecode::FLOAT;
  out.val = as_int (i);
  return (out.as_obj ());
}

static inline object
bigint2flt (object i, local_varobj<bigfloat>& out)
{
  const auto lp = as_bigint (i);
  out.local_init (lp->data, lp->len, F_ABS (lp->len));
#ifdef KP_ARCH_WIDE
  return (out.as_obj () | (i & SIGN_BIT));
#else
  return (out.as_obj ());
#endif
}

static object
downgrade_float (interpreter *interp, object obj)
{
  if (float_p (obj))
    {
      double i, val = as_float (obj);
      if (modf (val, &i) == 0)
        {
          int64_t uv = (int64_t)val;
          return (intobj (interp, uv));
        }
    }
  else
    {
      const auto fp = as_bigfloat (obj);
      int fl = F_ABS (fp->len);

      if (fl == fp->expo)
        {
          auto lp = as_bigint (alloc_bigint (interp, fl));
          memcpy (lp->data, fp->data, fl * sizeof (*fp->data));
          lp->len = fp->len;

#ifdef KP_ARCH_WIDE
          interp->alval |= obj & SIGN_BIT;
#endif
          return (interp->alval);
        }
    }

  return (obj);
}

#define MIXED_FN_NOCHECK()   ((void)0)
#define MIXED_FN_CHECK()   \
  if (ret.error_p ())   \
    return (exception ());   \
  else (void)0

#define MIXED_FN(name, ret_ty, check)   \
ret_ty name##_if (interpreter *interp, object x, object y)   \
{   \
  local_varobj<fltobj> tmp;   \
  ret_ty ret = name##_ff (interp, int2flt (x, tmp), y);   \
  check ();   \
  kp_return (downgrade_float (interp, deref (ret)));   \
}   \
\
object name##_iF (interpreter *interp, object x, object y)   \
{   \
  local_varobj<fltobj> tmp;   \
  kp_return (downgrade_float (interp,   \
                              name##_fF (interp, int2flt (x, tmp), y)));   \
}

MIXED_FN (add, object, MIXED_FN_NOCHECK)
MIXED_FN (sub, object, MIXED_FN_NOCHECK)
MIXED_FN (mul, object, MIXED_FN_NOCHECK)

#define MIXED_FN2(name)   \
MIXED_FN (name, result<object>, MIXED_FN_CHECK)   \
\
result<object> name##_fi (interpreter *interp, object x, object y)   \
{   \
  local_varobj<fltobj> tmp;   \
  object ret = KP_TRY (name##_ff (interp, x, int2flt (y, tmp)));   \
  kp_return (downgrade_float (interp, ret));   \
}   \
\
result<object> name##_Fi (interpreter *interp, object x, object y)   \
{   \
  local_varobj<fltobj> tmp;   \
  object ret = KP_TRY (name##_Ff (interp, x, int2flt (y, tmp)));   \
  kp_return (downgrade_float (interp, ret));   \
}

MIXED_FN2 (div)
MIXED_FN2 (mod)

#undef MIXED_FN
#undef MIXED_FN2

#define MIXED_FN(name, ret_ty, check)   \
ret_ty name##_If (interpreter *interp, object x, object y)   \
{   \
  int s2, c2;   \
  double v2 = get_dbl (y, s2, c2);   \
  if (c2 != FP_NORMAL)   \
    kp_return (y);   \
  \
  check (v2);   \
  local_varobj<bigfloat> tmp;   \
  fake_bigfloat ry (v2);   \
  \
  kp_return (downgrade_float (interp,   \
                              name##_FF (interp, bigint2flt (x, tmp),   \
                                         ry.get_obj (s2))));   \
}   \
\
object name##_IF (interpreter *interp, object x, object y)   \
{   \
  local_varobj<bigfloat> tmp;   \
  kp_return (downgrade_float (interp,   \
                              name##_FF (interp, bigint2flt (x, tmp), y)));   \
}   \

#undef MIXED_FN_NOCHECK
#define MIXED_FN_NOCHECK(x)   ((void)0)

MIXED_FN (add, object, MIXED_FN_NOCHECK)
MIXED_FN (sub, object, MIXED_FN_NOCHECK)
MIXED_FN (mul, object, MIXED_FN_NOCHECK)

#define MIXED_FN2(name)   \
MIXED_FN (name, result<object>, MIXED_FN_CHECK)   \
\
object name##_fI (interpreter *interp, object x, object y)   \
{   \
  int s1, c1;   \
  double v1 = get_dbl (x, s1, c1);   \
  if (c1 != FP_NORMAL)   \
    kp_return (x);   \
  \
  local_varobj<bigfloat> tmp;   \
  fake_bigfloat rx (v1);   \
  \
  kp_return (downgrade_float (interp,   \
                              name##_FF (interp, rx.get_obj (s1),   \
                                         bigint2flt (y, tmp))));   \
}   \
\
object name##_FI (interpreter *interp, object x, object y)   \
{   \
  local_varobj<bigfloat> tmp;   \
  kp_return (downgrade_float (interp,   \
                              name##_FF (interp, x, bigint2flt (y, tmp))));   \
}

#undef MIXED_FN_CHECK
#define MIXED_FN_CHECK(x)   \
  if ((x) == 0)   \
    return (interp->raise ("arith-error", "division by zero"));   \
  else (void)0

MIXED_FN2 (div)
MIXED_FN2 (mod)

// Division of integers that may produce floating point values.

result<object> div_ii (interpreter *interp, object x, object y)
{
  int v1 = as_int (x), v2 = as_int (y);
  if (v2 == 0)
    return (interp->raise ("arith-error", "division by zero"));
  kp_return (v1 % v2 == 0 ? fixint (v1 / v2) :
                            fltobj::make (interp, (double)v1 / v2));
}

object div_iI (interpreter *interp, object x, object y)
{
  local_varobj<fltobj> flt;
  local_varobj<bigfloat> bigf;
  return (div_fF (interp, int2flt (x, flt), bigint2flt (y, bigf)));
}

result<object> div_Ii (interpreter *interp, object x, object y)
{
  if (as_int (y) == 0)
    return (interp->raise ("arith-error", "division by zero"));

  local_varobj<fltobj> flt;
  local_varobj<bigfloat> bigf;
  return (downgrade_float (interp,
                           div_Ff (interp, bigint2flt (x, bigf),
                                   int2flt (y, flt))));
}

object div_II (interpreter *interp, object x, object y)
{
  local_varobj<bigfloat> lf1, lf2;
  return (downgrade_float (interp,
                           div_FF (interp, bigint2flt (x, lf1),
                                   bigint2flt (y, lf2))));
}

#undef MIXED_FN
#undef MIXED_FN2

int cmp_if (interpreter *interp, object x, object y)
{
  local_varobj<fltobj> fv;
  return (cmp_ff (interp, int2flt (x, fv), y));
}

int cmp_iF (interpreter *interp, object x, object y)
{
  local_varobj<fltobj> fv;
  return (cmp_fF (interp, int2flt (x, fv), y));
}

int cmp_If (interpreter *interp, object x, object y)
{
  local_varobj<bigfloat> fv;
  return (-cmp_fF (interp, y, bigint2flt (x, fv)));
}

int cmp_IF (interpreter *interp, object x, object y)
{
  local_varobj<bigfloat> fv;
  return (cmp_FF (interp, bigint2flt (x, fv), y));
}

bool eq_if (interpreter *, object x, object y)
{
  return ((double)as_int (x) == as_float (y));
}

bool eq_If (interpreter *, object x, object y)
{
  const auto lp = as_bigint (x);
  int xl = F_ABS (lp->len);

  if (xl > LNDBL_SIZE)
    return (false);

  double dbl = uitodbl (lp->data, xl, xl);
#ifdef KP_ARCH_WIDE
  if (x & SIGN_BIT)
#else
  if (lp->len < 0)
#endif
    dbl = -dbl;

  return (as_float (y) == dbl);
}

bool eq_IF (interpreter *, object x, object y)
{
  const auto lp = as_bigint (x);
  const auto fp = as_bigfloat (y);

#ifdef KP_ARCH_WIDE
  if ((x & SIGN_BIT) != (y & SIGN_BIT))
#else
  if ((lp->len ^ fp->len) < 0)
#endif
    return (false);

  return (F_ABS (fp->len) == fp->expo && fp->len == lp->len &&
          memcmp (fp->data, lp->data, F_ABS (lp->len) * sizeof (limb_t)) == 0);
}

object lsh_fi (interpreter *interp, object x, object y)
{
  int s1, c1, ival = as_int (y);
  double d1 = get_dbl (x, s1, c1);

  if (kp_unlikely (c1 != FP_NORMAL))
    kp_return (x);   // +/- Inf or NaN

  double ret = ldexp (d1, ival);
  if (kp_likely (!finf_p (ret)))
    kp_return (fltobj::make (interp, ret));

  fake_bigfloat tx (d1);
  return (lsh_Fi (interp, tx.as_obj (), y));
}

object lsh_Fi (interpreter *interp, object x, object y)
{
  int sign, ival = as_int (y);
  if (ival < 0)
    return (rsh_Fi (interp, x, -ival));

  const auto flt = get_bigfloat (x, sign);
  int len = F_ABS (flt->len);
  bigfloat *ret;

  if (ival % LIMB_BITS == 0)
    {
      ret = bigfloat::alloc_raw (len);
      memcpy (ret->data, flt->data, len * sizeof (limb_t));
      ret->expo = flt->expo + flt->expo / LIMB_BITS;
      ret->len = len;
    }
  else
    {
      ret = bigfloat::alloc_raw (len + 1);
      limb_t cy = uimul2exp (ret->data, flt->data,
                             len, ival % LIMB_BITS);

      ret->data[len] = cy;
      ret->len = len + (cy != 0);
      ret->expo = flt->expo + ival / LIMB_BITS + (cy != 0);
    }

  return (ret_F (interp, ret, sign));
}

object rsh_fi (interpreter *interp, object x, object y)
{
  int s1, c1, ival = as_int (y);
  double d1 = get_dbl (x, s1, c1);

  if (kp_unlikely (c1 != FP_NORMAL))
    kp_return (x);   // +/- Inf or NaN

  double ret = ldexp (d1, -ival);
  if (kp_likely (!finf_p (ret)))
    kp_return (fltobj::make (interp, ret));

  fake_bigfloat tx (d1);
  return (rsh_Fi (interp, tx.as_obj (), y));
}

object rsh_Fi (interpreter *interp, object x, object y)
{
  int sign, ival = as_int (y);
  if (ival < 0)
    return (rsh_Fi (interp, x, -ival));

  const auto flt = get_bigfloat (x, sign);
  int len = F_ABS (flt->len);
  bigfloat *ret;

  if (ival % LIMB_BITS == 0)
    {
      ret = bigfloat::alloc_raw (len);
      memcpy (ret->data, flt->data, len * sizeof (limb_t));
      ret->expo = flt->expo - flt->expo / LIMB_BITS;
      ret->len = len;
    }
  else
    {
      ret = bigfloat::alloc_raw (len + 1);
      limb_t cy = uidiv2exp (ret->data, flt->data,
                             len, LIMB_BITS - ival % LIMB_BITS);

      ret->data[len] = cy;
      ret->len = len + (cy != 0);
      ret->expo = flt->expo - ival / LIMB_BITS + (cy != 0) - 1;
    }

  return (ret_F (interp, ret, sign));
}

// External definitions.
object FLT_PINF;
object FLT_NINF;
object FLT_QNAN;
object FLT_ZERO;

static int
do_init_float (interpreter *)
{
  static fltobj SF_PINF;
  static fltobj SF_QNAN;
  static fltobj SF_ZERO;

#ifndef KP_ARCH_WIDE
  static fltobj SF_NINF;
  SF_NINF.vo_type = typecode::FLOAT;
#endif

  SF_PINF.vo_type = SF_QNAN.vo_type = SF_ZERO.vo_type = typecode::FLOAT;
  SF_PINF.val = std::numeric_limits<double>::infinity ();
  SF_QNAN.val = std::numeric_limits<double>::quiet_NaN ();
  SF_ZERO.val = 0.;

  // Now export the global objects.
#ifdef KP_ARCH_WIDE
  /* On 64-bit platforms, we use an additional bit to store the sign.
   * This helps us to avoid allocating a full value when negating it. */

#define DEF_FLOAT(ptr)   \
    ptrtype (ensure_mask ((ptr), TYPE_SHIFT - 1), typecode::FLOAT)

  FLT_PINF = DEF_FLOAT (&SF_PINF);
  FLT_NINF = FLT_PINF | SIGN_BIT;
  FLT_QNAN = DEF_FLOAT (&SF_QNAN);
  FLT_ZERO = DEF_FLOAT (&SF_ZERO);

#  undef DEF_FLOAT

#else
  SF_NINF.val = -SF_PINF.val;

  FLT_PINF = SF_PINF.as_obj ();
  FLT_NINF = SF_NINF.as_obj ();
  FLT_QNAN = SF_QNAN.as_obj ();
  FLT_ZERO = SF_ZERO.as_obj ();
#endif

  return (init_op::result_ok);
}

init_op init_float (do_init_float, "float");

KP_DECLS_END