/*  bignum.c -- bignum support                                */
/*  Copyright (c) 2009-2015 Alex Shinn.  All rights reserved. */
/*  BSD-style license: http://synthcode.com/license.txt       */

#include "chibi/sexp.h"

#if SEXP_USE_BIGNUMS

#define SEXP_INIT_BIGNUM_SIZE 2

static int digit_value (int c) {
  return (((c)<='9') ? ((c) - '0') : ((sexp_toupper(c) - 'A') + 10));
}

static int hex_digit (int n) {
  return ((n<=9) ? ('0' + n) : ('A' + n - 10));
}

sexp sexp_make_bignum (sexp ctx, sexp_uint_t len) {
  sexp_uint_t size = sexp_sizeof(bignum) + len*sizeof(sexp_uint_t);
  sexp res = sexp_alloc_tagged(ctx, size, SEXP_BIGNUM);
  if (!sexp_exceptionp(res)) {
    sexp_bignum_length(res) = len;
    sexp_bignum_sign(res) = 1;
  }
  return res;
}

sexp sexp_fixnum_to_bignum (sexp ctx, sexp a) {
  sexp res = sexp_make_bignum(ctx, 1);
  if (!sexp_exceptionp(res)) {
    sexp_bignum_data(res)[0] = sexp_unbox_fx_abs(a);
    sexp_bignum_sign(res) = sexp_fx_sign(a);
  }
  return res;
}

sexp sexp_make_integer_from_lsint (sexp ctx, sexp_lsint_t x) {
  sexp res;
  if (lsint_is_fixnum(x)) {
    res = sexp_make_fixnum(lsint_to_sint(x));
  } else if (sexp_lsint_fits_sint(x)) {
    res = sexp_make_bignum(ctx, 1);
    if (lsint_lt_0(x)) {
      sexp_bignum_sign(res) = -1;
      sexp_bignum_data(res)[0] = (sexp_uint_t)-lsint_to_sint(x);
    } else {
      sexp_bignum_sign(res) = 1;
      sexp_bignum_data(res)[0] = (sexp_uint_t)lsint_to_sint(x);
    }
  } else {
    res = sexp_make_bignum(ctx, 2);
    if (lsint_lt_0(x)) {
      sexp_bignum_sign(res) = -1;
      sexp_bignum_data(res)[0] = (sexp_uint_t)-lsint_to_sint(x);
      sexp_bignum_data(res)[1] = (sexp_uint_t)~lsint_to_sint_hi(x);
    } else {
      sexp_bignum_sign(res) = 1;
      sexp_bignum_data(res)[0] = (sexp_uint_t)lsint_to_sint(x);
      sexp_bignum_data(res)[1] = (sexp_uint_t)lsint_to_sint_hi(x);
    }
  }
  return res;
}

sexp sexp_make_unsigned_integer_from_luint (sexp ctx, sexp_luint_t x) {
  sexp res;
  if (luint_is_fixnum(x)) {
    res = sexp_make_fixnum(luint_to_uint(x));
  } else if (sexp_luint_fits_uint(x)) {
    res = sexp_make_bignum(ctx, 1);
    sexp_bignum_sign(res) = 1;
    sexp_bignum_data(res)[0] = luint_to_uint(x);
  } else {
    res = sexp_make_bignum(ctx, 2);
    sexp_bignum_sign(res) = 1;
    sexp_bignum_data(res)[0] = luint_to_uint(x);
    sexp_bignum_data(res)[1] = luint_to_uint_hi(x);
  }
  return res;
}

#if SEXP_USE_CUSTOM_LONG_LONGS
sexp sexp_make_integer(sexp ctx, long long x) {
  return sexp_make_integer_from_lsint(ctx, lsint_from_sint(x));
}
sexp sexp_make_unsigned_integer(sexp ctx, unsigned long long x) {
  return sexp_make_unsigned_integer_from_luint(ctx, luint_from_uint(x));
}
#else
sexp sexp_make_integer (sexp ctx, sexp_lsint_t x) {
  return sexp_make_integer_from_lsint(ctx, x);
}
sexp sexp_make_unsigned_integer (sexp ctx, sexp_luint_t x) {
  return sexp_make_unsigned_integer_from_luint(ctx, x);
}
#endif

#if !SEXP_64_BIT
long long sexp_bignum_to_sint(sexp x) {
  if (!sexp_bignump(x))
    return 0;
  if (sexp_bignum_length(x) > 1)
    return sexp_bignum_sign(x) * (
      (((long long)(sexp_bignum_data(x)[1]))<<(8*sizeof(sexp_bignum_data(x)[0]))) + sexp_bignum_data(x)[0]);
  return sexp_bignum_sign(x) * sexp_bignum_data(x)[0];
}

unsigned long long sexp_bignum_to_uint(sexp x) {
  if (!sexp_bignump(x))
    return 0;
  if (sexp_bignum_length(x) > 1)
    return (((unsigned long long)(sexp_bignum_data(x)[1]))<<(8*sizeof(sexp_bignum_data(x)[0]))) + sexp_bignum_data(x)[0];
  return sexp_bignum_data(x)[0];
}
#endif

#define double_trunc_10s_digit(f) (trunc((f)/10.0)*10.0)
#define double_10s_digit(f) ((f)-double_trunc_10s_digit(f))

#define double_16s_digit(f) fmod(f,16.0)

sexp sexp_double_to_bignum (sexp ctx, double f) {
  int sign;
  sexp_gc_var3(res, scale, tmp);
  sexp_gc_preserve3(ctx, res, scale, tmp);
  res = sexp_fixnum_to_bignum(ctx, SEXP_ZERO);
  scale = sexp_fixnum_to_bignum(ctx, SEXP_ONE);
  sign = (f < 0 ? -1 : 1);
  for (f=fabs(f); f >= 1.0; f=trunc(f/16)) {
    tmp = sexp_bignum_fxmul(ctx, NULL, scale, (sexp_uint_t)double_16s_digit(f), 0);
    res = sexp_bignum_add(ctx, res, res, tmp);
    scale = sexp_bignum_fxmul(ctx, NULL, scale, 16, 0);
  }
  sexp_bignum_sign(res) = sign;
  sexp_gc_release3(ctx);
  return res;
}

sexp sexp_copy_bignum (sexp ctx, sexp dst, sexp a, sexp_uint_t len0) {
  sexp_uint_t len = (len0 > 0) ? len0 : sexp_bignum_length(a), size;
  size = sexp_sizeof(bignum) + len*sizeof(sexp_uint_t);
  if (! dst || sexp_bignum_length(dst) < len) {
    dst = sexp_alloc_tagged(ctx, size, SEXP_BIGNUM);
    sexp_bignum_length(dst) = len;
  }
  if (sexp_bignum_length(a) < len)
    len = sexp_bignum_length(a);

  sexp_bignum_sign(dst) = sexp_bignum_sign(a);
  memset(sexp_bignum_data(dst), 0,
         sexp_bignum_length(dst)*sizeof(sexp_uint_t));
  memmove(sexp_bignum_data(dst), sexp_bignum_data(a),
          len*sizeof(sexp_uint_t));
  return dst;
}

int sexp_bignum_zerop (sexp a) {
  int i;
  sexp_uint_t *data = sexp_bignum_data(a);
  for (i=sexp_bignum_length(a)-1; i>=0; i--)
    if (data[i])
      return 0;
  return 1;
}

sexp_uint_t sexp_bignum_hi (sexp a) {
  sexp_uint_t i=sexp_bignum_length(a)-1;
  while ((i>0) && ! sexp_bignum_data(a)[i])
    i--;
  return i+1;
}

sexp_sint_t sexp_bignum_compare_abs (sexp a, sexp b) {
  int ai=sexp_bignum_hi(a), bi=sexp_bignum_hi(b);
  sexp_uint_t *adata=sexp_bignum_data(a), *bdata=sexp_bignum_data(b);
  if (ai != bi)
    return ai - bi;
  for (--ai; ai >= 0; ai--) {
    if (adata[ai] > bdata[ai])
      return 1;
    else if (adata[ai] < bdata[ai])
      return -1;
  }
  return 0;
}

sexp_sint_t sexp_bignum_compare (sexp a, sexp b) {
  if (sexp_bignum_sign(a) != sexp_bignum_sign(b))
    return sexp_bignum_sign(a);
  sexp_sint_t cmp = sexp_bignum_compare_abs(a, b);
  return sexp_bignum_sign(a) < 0 ? -cmp : cmp;
}

sexp sexp_bignum_normalize (sexp a) {
  sexp_uint_t *data;
  if ((! sexp_bignump(a)) || (sexp_bignum_hi(a)>1))
    return a;
  data = sexp_bignum_data(a);
  if ((data[0] > SEXP_MAX_FIXNUM)
      && ! ((sexp_bignum_sign(a) == -1) && (data[0] == SEXP_MAX_FIXNUM+1)))
    return a;
  return sexp_make_fixnum((sexp_sint_t)data[0] * sexp_bignum_sign(a));
}

double sexp_bignum_to_double (sexp a) {
  double res = 0;
  sexp_sint_t i;
  sexp_uint_t *data=sexp_bignum_data(a);
  for (i=sexp_bignum_hi(a)-1; i>=0; i--)
    res = res * ((double)SEXP_UINT_T_MAX+1) + data[i];
  return res * sexp_bignum_sign(a);
}

sexp sexp_bignum_fxadd (sexp ctx, sexp a, sexp_uint_t b) {
  sexp_uint_t len=sexp_bignum_hi(a), *data=sexp_bignum_data(a),
    carry=b, i=0, n;
  do { n = data[i];
       data[i] += carry;
       carry = (n > (SEXP_UINT_T_MAX - carry));
  } while (++i<len && carry);
  if (carry) {
    a = sexp_copy_bignum(ctx, NULL, a, len+1);
    sexp_bignum_data(a)[len] = 1;
  }
  return a;
}

sexp sexp_bignum_fxsub (sexp ctx, sexp a, sexp_uint_t b) {
  sexp_uint_t len=sexp_bignum_hi(a), *data=sexp_bignum_data(a), borrow, i=0, n;
  if ((len == 1) && (b > data[0])) {
    data[0] = b - data[0];
    sexp_bignum_sign(a) = -sexp_bignum_sign(a);
  } else {
    for (borrow=b; borrow; i++) {
      n = data[i];
      data[i] -= borrow;
      borrow = (n < borrow);
    }
  }
  return a;
}

sexp sexp_bignum_fxmul (sexp ctx, sexp d, sexp a, sexp_uint_t b, int offset) {
  sexp_uint_t len=sexp_bignum_length(a), *data, *adata=sexp_bignum_data(a),
    carry=0, i;
  sexp_luint_t n;
  sexp_gc_var1(tmp);
  sexp_gc_preserve1(ctx, tmp);
  if ((! d) || (sexp_bignum_length(d) < len+offset))
    d = sexp_make_bignum(ctx, len+offset);
  tmp = d;
  data = sexp_bignum_data(d);
  for (i=0; i<len; i++) {
    n = luint_add(luint_mul_uint(luint_from_uint(adata[i]), b), luint_from_uint(carry));
    data[i+offset] = luint_to_uint(n);
    carry = luint_to_uint(luint_shr(n, (sizeof(sexp_uint_t)*8)));
  }
  if (carry) {
    if (sexp_bignum_length(d) <= len+offset)
      d = sexp_copy_bignum(ctx, NULL, d, len+offset+1);
    sexp_bignum_data(d)[len+offset] = carry;
  }
  sexp_gc_release1(ctx);
  return d;
}

sexp_uint_t sexp_bignum_fxdiv (sexp ctx, sexp a, sexp_uint_t b, int offset) {
  sexp_uint_t len=sexp_bignum_hi(a), *data=sexp_bignum_data(a), q, r=0;
  int i;
  sexp_luint_t n = luint_from_uint(0);
  for (i=len-1; i>=offset; i--) {
    n = luint_add(luint_shl(n,  sizeof(sexp_uint_t)*8), luint_from_uint(data[i]));
    q = luint_to_uint(luint_div_uint(n, b));
    r = luint_to_uint(luint_sub(n, luint_mul_uint(luint_from_uint(q),  b)));
    data[i] = q;
    n = luint_from_uint(r);
  }
  return r;
}

sexp sexp_bignum_fxrem (sexp ctx, sexp a, sexp_sint_t b) {
  sexp_uint_t len=sexp_bignum_hi(a), *data=sexp_bignum_data(a), q, b0;
  int i;
  sexp_luint_t n = luint_from_uint(0);
  if (b > 0) {
    q = b - 1;
    if ((b & q) == 0)
      return sexp_make_fixnum(sexp_bignum_sign(a) * (data[0] & q));
  }
  b0 = (b >= 0) ? b : -b;
  if (b0 == 0) {
    return sexp_xtype_exception(ctx, NULL, "divide by zero", a);
  }
  for (i=len-1; i>=0; i--) {
    n = luint_add(luint_shl(n, sizeof(sexp_uint_t)*8), luint_from_uint(data[i]));
    q = luint_to_uint(luint_div_uint(n, b0));
    n = luint_sub(n, luint_mul_uint(luint_from_uint(q), b0));
  }
  return sexp_make_fixnum(sexp_bignum_sign(a) * (sexp_sint_t)luint_to_uint(n));
}

sexp sexp_read_bignum (sexp ctx, sexp in, sexp_uint_t init,
                       signed char sign, sexp_uint_t base) {
  int c, digit;
  sexp_gc_var3(res, tmp, imag);
  sexp_gc_preserve3(ctx, res, tmp, imag);
  res = sexp_make_bignum(ctx, SEXP_INIT_BIGNUM_SIZE);
  sexp_bignum_sign(res) = sign;
  sexp_bignum_data(res)[0] = init;
  for (c=sexp_read_char(ctx, in); sexp_isxdigit(c); c=sexp_read_char(ctx, in)) {
    digit = digit_value(c);
    if ((digit < 0) || (digit >= (int)base))
      break;
    res = sexp_bignum_fxmul(ctx, res, res, base, 0);
    res = sexp_bignum_fxadd(ctx, res, digit);
  }
  if (c=='.' || c=='e' || c=='E') {
    if (base != 10) {
      res = sexp_read_error(ctx, "found non-base 10 float", SEXP_NULL, in);
    } else if (c=='.') {
      res = sexp_read_float_tail(ctx, in, sexp_bignum_to_double(res), (sign==-1));
    } else {
      tmp = sexp_read_number(ctx, in, base, 0);
#if SEXP_USE_COMPLEX
      if (sexp_complexp(tmp)) {
        imag = sexp_complex_imag(tmp);
        tmp = sexp_complex_real(tmp);
      } else {
        imag = SEXP_ZERO;
      }
#endif
      if (sexp_exceptionp(tmp)) {
        res = tmp;
      } else if (sexp_fixnump(tmp) && labs(sexp_unbox_fixnum(tmp)) < 100*1024*1024) {
        tmp = sexp_expt(ctx, SEXP_TEN, tmp);
        res = sexp_mul(ctx, res, tmp);
      } else {
        tmp = sexp_exact_to_inexact(ctx, NULL, 2, tmp);
        tmp = sexp_expt(ctx, SEXP_TEN, tmp);
        res = sexp_mul(ctx, res, tmp);
      }
#if SEXP_USE_COMPLEX
      if (imag != SEXP_ZERO && !sexp_exceptionp(res))
        res = sexp_make_complex(ctx, res, imag);
#endif
    }
#if SEXP_USE_RATIOS
  } else if (c=='/') {
    res = sexp_bignum_normalize(res);
    res = sexp_make_ratio(ctx, res, SEXP_ONE);
    sexp_ratio_denominator(res) = sexp_read_number(ctx, in, 10, 0);
    res = sexp_ratio_normalize(ctx, res, in);
#endif
#if SEXP_USE_COMPLEX
  } else if (c=='i' || c=='i' || c=='+' || c=='-') {
    sexp_push_char(ctx, c, in);
    res = sexp_bignum_normalize(res);
    res = sexp_read_complex_tail(ctx, in, res);
#endif
  } else if ((c!=EOF) && ! sexp_is_separator(c)) {
    res = sexp_read_error(ctx, "invalid numeric syntax",
                          sexp_make_character(c), in);
  } else {
    sexp_push_char(ctx, c, in);
  }
  sexp_gc_release3(ctx);
  return sexp_bignum_normalize(res);
}

static int log2i(int v) {
  int i;
  for (i = 0; i < sizeof(v)*8; i++)
    if ((1<<(i+1)) > v)
      break;
  return i;
}

sexp sexp_write_bignum (sexp ctx, sexp a, sexp out, sexp_uint_t base) {
  int i, str_len, lg_base = log2i(base);
  char *data;
  sexp_gc_var2(b, str);
  sexp_gc_preserve2(ctx, b, str);
  b = sexp_copy_bignum(ctx, NULL, a, 0);
  sexp_bignum_sign(b) = 1;
  if (lg_base < 1) {
    return sexp_xtype_exception(ctx, NULL, "number base too small", a);
  }
  i = str_len = (sexp_bignum_length(b)*sizeof(sexp_uint_t)*8 + lg_base - 1)
    / lg_base + 1;
  str = sexp_make_string(ctx, sexp_make_fixnum(str_len),
                         sexp_make_character(' '));
  data = sexp_string_data(str);
  while (! sexp_bignum_zerop(b))
    data[--i] = hex_digit(sexp_bignum_fxdiv(ctx, b, base, 0));
  if (i == str_len)
    data[--i] = '0';
  else if (sexp_bignum_sign(a) == -1)
    data[--i] = '-';
  sexp_write_string(ctx, data + i, out);
  sexp_gc_release2(ctx);
  return SEXP_VOID;
}

/****************** bignum arithmetic *************************/

sexp sexp_bignum_add_fixnum (sexp ctx, sexp a, sexp b) {
  sexp_gc_var1(c);
  sexp_gc_preserve1(ctx, c);
  c = sexp_copy_bignum(ctx, NULL, a, 0);
  if (sexp_bignum_sign(c) == sexp_fx_sign(b))
    c = sexp_bignum_fxadd(ctx, c, sexp_unbox_fx_abs(b));
  else
    c = sexp_bignum_fxsub(ctx, c, sexp_unbox_fx_abs(b));
  sexp_gc_release1(ctx);
  return c;
}

sexp sexp_bignum_sub_digits (sexp ctx, sexp dst, sexp a, sexp b) {
  sexp_uint_t alen=sexp_bignum_hi(a), blen=sexp_bignum_hi(b),
    borrow=0, i, *adata, *bdata, *cdata;
  sexp_gc_var1(c);
  if ((alen < blen) || ((alen == blen) && (sexp_bignum_compare_abs(a, b) < 0)))
    return sexp_bignum_sub_digits(ctx, dst, b, a);
  sexp_gc_preserve1(ctx, c);
  c = ((dst && sexp_bignum_hi(dst) >= alen)
       ? dst : sexp_copy_bignum(ctx, NULL, a, 0));
  adata = sexp_bignum_data(a);
  bdata = sexp_bignum_data(b);
  cdata = sexp_bignum_data(c);
  for (i=0; i<blen; i++) {
    if (adata[i] > bdata[i] || (adata[i] == bdata[i] && !borrow)) {
      cdata[i] = adata[i] - bdata[i] - borrow;
      borrow = 0;
    } else {
      cdata[i] = (SEXP_UINT_T_MAX - bdata[i]);
      cdata[i] += 1;
      cdata[i] -= borrow;
      cdata[i] += adata[i];
      borrow = 1;
    }
  }
  for ( ; borrow && (i<alen); i++) {
    borrow = (cdata[i] == 0 ? 1 : 0);
    cdata[i]--;
  }
  sexp_gc_release1(ctx);
  return c;
}

sexp sexp_bignum_add_digits (sexp ctx, sexp dst, sexp a, sexp b) {
  sexp_uint_t alen=sexp_bignum_hi(a), blen=sexp_bignum_hi(b),
    carry=0, i, old_a, p_sum, *adata, *bdata, *cdata;
  sexp_gc_var1(c);
  if (alen < blen) return sexp_bignum_add_digits(ctx, dst, b, a);
  sexp_gc_preserve1(ctx, c);
  c = ((dst && sexp_bignum_hi(dst) >= alen)
       ? dst : sexp_copy_bignum(ctx, NULL, a, 0));
  adata = sexp_bignum_data(a);
  bdata = sexp_bignum_data(b);
  cdata = sexp_bignum_data(c);
  for (i=0; i<blen; i++) {
    old_a = adata[i]; /* adata may alias cdata */
    p_sum = adata[i] + bdata[i];
    cdata[i] = p_sum + carry;
    carry = (old_a > (SEXP_UINT_T_MAX - bdata[i]) ? 1 : 0)
          + (p_sum > (SEXP_UINT_T_MAX - carry) ? 1 : 0);
  }
  for ( ; carry && (i<alen); i++) {
    carry = (cdata[i] == SEXP_UINT_T_MAX ? 1 : 0);
    cdata[i]++;
  }
  if (carry) {
    c = sexp_copy_bignum(ctx, NULL, c, alen+1);
    sexp_bignum_data(c)[alen] = 1;
  }
  sexp_gc_release1(ctx);
  return c;
}

sexp sexp_bignum_add (sexp ctx, sexp dst, sexp a, sexp b) {
  sexp res;
  if (sexp_bignum_sign(a) == sexp_bignum_sign(b)) {
    res = sexp_bignum_add_digits(ctx, dst, a, b);
    sexp_bignum_sign(res) = sexp_bignum_sign(a);
  } else {
    res = sexp_bignum_sub_digits(ctx, dst, a, b);
    sexp_bignum_sign(res)
      = sexp_bignum_sign(sexp_bignum_compare_abs(a, b) >= 0 ? a : b);
  }
  return res;
}

sexp sexp_bignum_sub (sexp ctx, sexp dst, sexp a, sexp b) {
  sexp res;
  int sign;
  if (sexp_bignum_sign(a) == sexp_bignum_sign(b)) {
    sign = (sexp_bignum_compare_abs(a, b) >= 0 ? sexp_bignum_sign(a)
            : -sexp_bignum_sign(a));
    res = sexp_bignum_sub_digits(ctx, dst, a, b);
  } else {
    sign = sexp_bignum_sign(a);
    res = sexp_bignum_add_digits(ctx, dst, a, b);
  }
  sexp_bignum_sign(res) = sign;
  return res;
}

static void sexp_bignum_split (sexp ctx, sexp a, sexp_uint_t k, sexp* lo, sexp* hi) {
  sexp_uint_t alen=sexp_bignum_hi(a), i, *adata=sexp_bignum_data(a),
    *lodata, *hidata;
  *lo = sexp_make_bignum(ctx, k);        /* must be gc protected by caller */
  *hi = sexp_make_bignum(ctx, alen-k+1);
  lodata = sexp_bignum_data(*lo);
  hidata = sexp_bignum_data(*hi);
  for (i=0; i<k; i++)                    /* split into a[0..k-1], a[k..] */
    lodata[i] = adata[i];
  for (i=k; i<alen; i++)
    hidata[i-k] = adata[i];
}

static sexp sexp_bignum_shift (sexp ctx, sexp a, sexp_uint_t k) {
  sexp res;
  sexp_uint_t alen = sexp_bignum_hi(a), i;
  res = sexp_make_bignum(ctx, alen + k + 1);
  for (i=0; i<alen; i++)
    sexp_bignum_data(res)[i+k] = sexp_bignum_data(a)[i];
  return res;
}

sexp sexp_bignum_mul (sexp ctx, sexp dst, sexp a, sexp b) {
  sexp_uint_t alen=sexp_bignum_hi(a), blen=sexp_bignum_hi(b), k,
    *bdata=sexp_bignum_data(b);
  sexp_gc_var7(a0, a1, b0, b1, z0, z1, z2);
  if (alen < blen) return sexp_bignum_mul(ctx, dst, b, a);
  if (blen == 1) {
    z1 = sexp_bignum_fxmul(ctx, dst, a, bdata[0], 0);
  } else {
    /* karatsuba:                             */
    /*   ab = (a1B^k + a0) * (b1B^k + b0)     */
    /*      = z2B^2k + z1B^k + z0             */
    /* where:                                 */
    /*   z2 = a1b1                            */
    /*   z1 = a1b0 + a0b1                     */
    /*   z0 = a0b0                            */
    /* then optimize further:                 */
    /*   z1 = (a1 + a0)(b1 + b0) - z2 - z0    */
    sexp_gc_preserve7(ctx, a0, a1, b0, b1, z0, z1, z2);
    k = blen / 2;
    sexp_bignum_split(ctx, a, k, &a0, &a1);
    sexp_bignum_split(ctx, b, k, &b0, &b1);
    z0 = sexp_bignum_add(ctx, NULL, a1, a0);  /* temp */
    z1 = sexp_bignum_add(ctx, NULL, b1, b0);
    z1 = sexp_bignum_mul(ctx, NULL, z1, z0);  /* 1 */
    z0 = sexp_bignum_mul(ctx, NULL, a0, b0);  /* 2 */
    z2 = sexp_bignum_mul(ctx, NULL, a1, b1);  /* 3 */
    z1 = sexp_bignum_sub(ctx, NULL, z1, z0);
    z1 = sexp_bignum_sub(ctx, NULL, z1, z2);
    z2 = sexp_bignum_shift(ctx, z2, 2*k);
    z1 = sexp_bignum_shift(ctx, z1, k);
    z1 = sexp_bignum_add(ctx, z1, z1, z0);
    z1 = sexp_bignum_add(ctx, z1, z1, z2);
    sexp_gc_release7(ctx);
  }
  sexp_bignum_sign(z1) = sexp_bignum_sign(a) * sexp_bignum_sign(b);
  return z1;
}

sexp sexp_bignum_quot_rem (sexp ctx, sexp *rem, sexp a, sexp b) {
  sexp_luint_t d, dn, dd;
  sexp_uint_t dlo, dhi;
  sexp_sint_t alen, blen=sexp_bignum_hi(b), sign=1, off=0;
  sexp_gc_var5(q, x, y, a1, b1);
  if (blen == 1 && sexp_bignum_data(b)[0] == 0)
    return sexp_xtype_exception(ctx, NULL, "divide by zero", a);
  sexp_gc_preserve5(ctx, q, x, y, a1, b1);
  a1 = sexp_copy_bignum(ctx, NULL, a, 0);
  sexp_bignum_sign(a1) = 1;
  /* fast path for single bigit divisor */
  if (blen == 1) {
    b1 = sexp_make_bignum(ctx, 1);
    sexp_bignum_data(b1)[0] = sexp_bignum_fxdiv(ctx, a1, sexp_bignum_data(b)[0], 0);
    *rem = sexp_bignum_normalize(b1);
    if (sexp_bignum_sign(a) * sexp_bignum_sign(b) < 0) {
      sexp_negate_exact(a1);
    }
    if (sexp_bignum_sign(a) < 0) {
      sexp_negate_exact(*rem);
    }
    sexp_gc_release5(ctx);
    return a1;
  }
  /* general case */
  b1 = sexp_copy_bignum(ctx, NULL, b, 0);
  sexp_bignum_sign(b1) = 1;
  q = SEXP_ZERO;
  x = sexp_make_bignum(ctx, sexp_bignum_length(a));
  while (sexp_bignum_compare_abs(a1, b1) >= 0) { /* a1, b1 at least 2 bigits */
    /* guess divisor x */
    alen = sexp_bignum_hi(a1);
    sexp_bignum_data(x)[off] = 0;
    if (off > 0) sexp_bignum_data(x)[off-1] = 0;
    off = alen - blen + 1;
    dn = luint_add_uint(luint_shl(luint_from_uint(sexp_bignum_data(a1)[alen-1])
           , (sizeof(sexp_uint_t)*8))
          , sexp_bignum_data(a1)[alen-2]);
    dd = luint_add_uint(luint_shl(luint_from_uint(sexp_bignum_data(b1)[blen-1])
           , (sizeof(sexp_uint_t)*8))
          , sexp_bignum_data(b1)[blen-2]);
    if (alen > 2 && blen > 2 &&
        luint_lt(luint_from_uint(sexp_bignum_data(a1)[alen-1]), luint_shl(luint_from_uint(1), (sizeof(sexp_uint_t)*4))) &&
        luint_lt(luint_from_uint(sexp_bignum_data(b1)[blen-1]), luint_shl(luint_from_uint(1), (sizeof(sexp_uint_t)*4)))) {
      dn = luint_add_uint(luint_shl(dn, (sizeof(sexp_uint_t)*4))
        , (sexp_bignum_data(a1)[alen-3] >> (sizeof(sexp_uint_t)*4)));
      dd = luint_add_uint(luint_shl(dd, (sizeof(sexp_uint_t)*4))
        , (sexp_bignum_data(b1)[blen-3] >> (sizeof(sexp_uint_t)*4)));
    }
    d = luint_div(dn, dd);
    if (luint_eq(d, luint_from_uint(0))) {
      dn = luint_add_uint(luint_shl(luint_from_uint(sexp_bignum_data(a1)[alen-1])
             , (sizeof(sexp_uint_t)*8))
            , sexp_bignum_data(a1)[alen-2]);
      dd = luint_from_uint(sexp_bignum_data(b1)[blen-1]);
      if (luint_lt(luint_from_uint(sexp_bignum_data(a1)[alen-1]), (luint_shl(luint_from_uint(1), (sizeof(sexp_uint_t)*4)))) &&
          luint_lt(luint_from_uint(sexp_bignum_data(b1)[blen-1]), (luint_shl(luint_from_uint(1), (sizeof(sexp_uint_t)*4))))) {
        dn = luint_add_uint(luint_shl(dn, (sizeof(sexp_uint_t)*4))
          , (sexp_bignum_data(a1)[alen-3] >> (sizeof(sexp_uint_t)*4)));
        dd = luint_add_uint(luint_shl(dd, (sizeof(sexp_uint_t)*4))
          , (sexp_bignum_data(b1)[blen-2] >> (sizeof(sexp_uint_t)*4)));
      }
      d = luint_div(dn, dd);
      off--;
    }
    dhi = luint_to_uint(luint_shr(d, (sizeof(sexp_uint_t)*8)));
    dlo = luint_to_uint(luint_and(d, luint_sub(luint_shl(luint_from_uint(1), (sizeof(sexp_uint_t)*8)), luint_from_uint(1))));
    sexp_bignum_data(x)[off] = dhi;
    if (off > 0) sexp_bignum_data(x)[off-1] = dlo;
    /* update quotient q and remainder a1 estimates */
    y = sexp_bignum_mul(ctx, NULL, b1, x);
    if (sign < 0) {
      a1 = sexp_bignum_add(ctx, NULL, a1, y);
      q = sexp_sub(ctx, q, x);
    } else {
      a1 = sexp_bignum_sub(ctx, NULL, a1, y);
      q = sexp_add(ctx, q, x);
    }
    /* flip the sign if we overshot in our estimate */
    if (sexp_bignum_sign(a1) != sign) {
      sexp_bignum_sign(a1) = (char)(-sign);
      sign *= -1;
    }
  }
  /* adjust signs */
  a1 = sexp_bignum_normalize(a1);
  if (sign < 0 && a1 != SEXP_ZERO) {
    q = sexp_sub(ctx, q, SEXP_ONE);
    a1 = sexp_add(ctx, a1, b1);
  }
  *rem = a1;
  if (sexp_bignum_sign(a) * sexp_bignum_sign(b) < 0) {
    sexp_negate_exact(q);
  }
  if (sexp_bignum_sign(a) < 0) {
    sexp_negate_exact(*rem);
  }
  sexp_gc_release5(ctx);
  return q;
}

sexp sexp_bignum_quotient (sexp ctx, sexp a, sexp b) {
  sexp res;
  sexp_gc_var1(rem);
  sexp_gc_preserve1(ctx, rem);
  res = sexp_bignum_quot_rem(ctx, &rem, a, b);
  sexp_gc_release1(ctx);
  return res;
}

sexp sexp_bignum_remainder (sexp ctx, sexp a, sexp b) {
  sexp_gc_var1(rem);
  sexp_gc_preserve1(ctx, rem);
  sexp_bignum_quot_rem(ctx, &rem, a, b); /* discard quotient */
  sexp_gc_release1(ctx);
  return rem;
}

sexp sexp_bignum_expt (sexp ctx, sexp a, sexp b) {
  sexp_sint_t e = sexp_unbox_fixnum(b);
  sexp_sint_t abs_e;
  if (e < 0)
    abs_e = -e;
  else
    abs_e = e;
  sexp_gc_var2(res, acc);
  sexp_gc_preserve2(ctx, res, acc);
  res = sexp_fixnum_to_bignum(ctx, SEXP_ONE);
  acc = sexp_copy_bignum(ctx, NULL, a, 0);
  for (; abs_e; abs_e>>=1, acc=sexp_bignum_mul(ctx, NULL, acc, acc))
    if (abs_e & 1)
      res = sexp_bignum_mul(ctx, NULL, res, acc);
  if (e < 0)
    res = sexp_div(ctx, sexp_fixnum_to_bignum(ctx, SEXP_ONE), res);
  sexp_gc_release2(ctx);
  return sexp_bignum_normalize(res);
}

#if SEXP_USE_MATH

/*
 * a = x * 2^2n, with 0.1 <= x < 1.0 (base 2)  =>  sqrt(a) ~ 2^n
 */
sexp sexp_bignum_sqrt_estimate (sexp ctx, sexp a) {
  sexp_uint_t alen=sexp_bignum_hi(a), adata_hi;
  int nbits, i;
  sexp_gc_var1(res);

  adata_hi = sexp_bignum_data(a)[alen - 1];
  for (i = sizeof(sexp_uint_t)*8-1; i > 0; i--)
    if (adata_hi & (1ul << i))
      break;
  nbits = sizeof(sexp_uint_t) * 8 * (alen - 1) + i + 1;

  sexp_gc_preserve1(ctx, res);
  res = sexp_fixnum_to_bignum(ctx, SEXP_TWO);
  res = sexp_bignum_expt(ctx, res, sexp_make_fixnum(nbits / 2));
  sexp_gc_release1(ctx);
  return res;
}

#define SEXP_MAX_ACCURATE_FLONUM_SQRT 1073741824.0 /* 2^30 */

/* Babylonian method */
sexp sexp_bignum_sqrt (sexp ctx, sexp a, sexp* rem_out) {
  sexp_gc_var4(res, rem, tmp, tmpa);
  if (! sexp_bignump(a)) return sexp_type_exception(ctx, NULL, SEXP_BIGNUM, a);
  sexp_gc_preserve4(ctx, res, rem, tmp, tmpa);
  /* initial estimate via flonum, ignoring signs */
  if (sexp_exact_negativep(a)) {
    tmpa = sexp_copy_bignum(ctx, NULL, a, 0);
    a = tmpa;
    sexp_negate(a);
  }
  res = sexp_make_flonum(ctx, sexp_bignum_to_double(a));
  res = sexp_inexact_sqrt(ctx, NULL, 1, res);
  if (sexp_flonump(res) &&
      sexp_flonum_value(res) > SEXP_MAX_ACCURATE_FLONUM_SQRT) {
    if (isinf(sexp_flonum_value(res)))
      res = sexp_bignum_sqrt_estimate(ctx, a);
    else
      res = sexp_double_to_bignum(ctx, sexp_flonum_value(res));
  loop:                           /* until 0 <= a - res*res < 2*res + 1 */
    rem = sexp_mul(ctx, res, res);
    tmp = rem = sexp_sub(ctx, a, rem);
    if (!sexp_positivep(tmp)) goto adjust;
    tmp = sexp_sub(ctx, tmp, SEXP_ONE);
    tmp = sexp_quotient(ctx, tmp, SEXP_TWO);
    tmp = sexp_compare(ctx, tmp, res);
    if (sexp_positivep(tmp)) {
    adjust:
      tmp = sexp_quotient(ctx, a, res);
      res = sexp_add(ctx, res, tmp);
      res = sexp_quotient(ctx, res, SEXP_TWO);
      goto loop;
    }
  } else {
    if (sexp_flonump(res))
      res = sexp_bignum_normalize(sexp_double_to_bignum(ctx, sexp_flonum_value(res)));
    tmp = sexp_mul(ctx, res, res);
    rem = sexp_sub(ctx, a, tmp);
  }
  *rem_out = sexp_bignum_normalize(rem);
  sexp_gc_release4(ctx);
  return sexp_bignum_normalize(res);
}

#endif  /* SEXP_USE_MATH */

/************************ ratios ******************************/

#if SEXP_USE_RATIOS

double sexp_ratio_to_double (sexp ctx, sexp rat) {
  sexp_gc_var1(quot);
  sexp num = sexp_ratio_numerator(rat), den = sexp_ratio_denominator(rat);
  double res = (sexp_bignump(num) ? sexp_bignum_to_double(num)
          : sexp_fixnum_to_double(num))
    / (sexp_bignump(den) ? sexp_bignum_to_double(den)
       : sexp_fixnum_to_double(den));
  if (!isfinite(res)) {
    sexp_gc_preserve1(ctx, quot);
    if (sexp_unbox_fixnum(sexp_compare(ctx, sexp_ratio_numerator(rat),  sexp_ratio_denominator(rat))) < 0) {
      quot = sexp_quotient(ctx, sexp_ratio_denominator(rat),  sexp_ratio_numerator(rat));
      res = 1 / sexp_to_double(ctx, quot);
    } else {
      quot = sexp_quotient(ctx, sexp_ratio_numerator(rat),  sexp_ratio_denominator(rat));
      res = sexp_to_double(ctx, quot);
    }
    sexp_gc_release1(ctx);
  }
  return res;
}

sexp sexp_double_to_ratio (sexp ctx, double f) {
  int sign, i;
  sexp_gc_var3(res, whole, scale);
  if (f == trunc(f))
    return sexp_bignum_normalize(sexp_double_to_bignum(ctx, f));
  sexp_gc_preserve3(ctx, res, whole, scale);
  whole = sexp_double_to_bignum(ctx, trunc(f));
  res = sexp_fixnum_to_bignum(ctx, SEXP_ZERO);
  scale = SEXP_ONE;
  sign = (f < 0 ? -1 : 1);
  for (i=0, f=fabs(f-trunc(f)); f != trunc(f) && i < 15; i++) {
    res = sexp_bignum_fxmul(ctx, NULL, res, 10, 0);
    f = f * 10;
    res = sexp_bignum_fxadd(ctx, res, (sexp_uint_t)double_10s_digit(f));
    f = f - trunc(f);
    scale = sexp_mul(ctx, scale, SEXP_TEN);
  }
  sexp_bignum_sign(res) = sign;
  res = sexp_bignum_normalize(res);
  scale = sexp_bignum_normalize(scale);
  res = sexp_make_ratio(ctx, res, scale);
  res = sexp_ratio_normalize(ctx, res, SEXP_FALSE);
  res = sexp_add(ctx, res, whole);
  sexp_gc_release3(ctx);
  return res;
}

/*
 * For conversion that does not introduce round-off error,
 * no matter what FLT_RADIX is.
 */
sexp sexp_double_to_ratio_2 (sexp ctx, double f) {
  int sign,i;
  sexp_gc_var3(res, whole, scale);
  if (f == trunc(f))
    return sexp_bignum_normalize(sexp_double_to_bignum(ctx, f));
  sexp_gc_preserve3(ctx, res, whole, scale);
  whole = sexp_double_to_bignum(ctx, trunc(f));
  res = sexp_fixnum_to_bignum(ctx, SEXP_ZERO);
  scale = SEXP_ONE;
  sign = (f < 0 ? -1 : 1);
  f = fabs(f-trunc(f));
  while(f) {
    res = sexp_bignum_fxmul(ctx, NULL, res, FLT_RADIX, 0);
    scale = sexp_mul(ctx, scale, sexp_make_fixnum(FLT_RADIX));
    f *= FLT_RADIX;
    i = trunc(f);
    if (i) {
      f -= i;
      res = sexp_bignum_fxadd(ctx, res, i);
    }
  }
  sexp_bignum_sign(res) = sign;
  res = sexp_bignum_normalize(res);
  scale = sexp_bignum_normalize(scale);
  res = sexp_make_ratio(ctx, res, scale);
  res = sexp_ratio_normalize(ctx, res, SEXP_FALSE);
  res = sexp_add(ctx, res, whole);
  sexp_gc_release3(ctx);
  return res;
}

sexp sexp_ratio_add (sexp ctx, sexp a, sexp b) {
  sexp_gc_var3(res, num, den);
  sexp_gc_preserve3(ctx, res, num, den);
  num = sexp_mul(ctx, sexp_ratio_numerator(a), sexp_ratio_denominator(b));
  den = sexp_mul(ctx, sexp_ratio_numerator(b), sexp_ratio_denominator(a));
  num = sexp_add(ctx, num, den);
  den = sexp_mul(ctx, sexp_ratio_denominator(a), sexp_ratio_denominator(b));
  res = sexp_make_ratio(ctx, num, den);
  res = sexp_ratio_normalize(ctx, res, SEXP_FALSE);
  sexp_gc_release3(ctx);
  return res;
}

sexp sexp_ratio_mul (sexp ctx, sexp a, sexp b) {
  sexp_gc_var3(res, num, den);
  sexp_gc_preserve3(ctx, res, num, den);
  num = sexp_mul(ctx, sexp_ratio_numerator(a), sexp_ratio_numerator(b));
  den = sexp_mul(ctx, sexp_ratio_denominator(a), sexp_ratio_denominator(b));
  res = sexp_make_ratio(ctx, num, den);
  res = sexp_ratio_normalize(ctx, res, SEXP_FALSE);
  sexp_gc_release3(ctx);
  return res;
}

sexp sexp_ratio_div (sexp ctx, sexp a, sexp b) {
  sexp_gc_var3(res, num, den);
  sexp_gc_preserve3(ctx, res, num, den);
  num = sexp_mul(ctx, sexp_ratio_numerator(a), sexp_ratio_denominator(b));
  den = sexp_mul(ctx, sexp_ratio_denominator(a), sexp_ratio_numerator(b));
  res = sexp_make_ratio(ctx, num, den);
  res = sexp_ratio_normalize(ctx, res, SEXP_FALSE);
  sexp_gc_release3(ctx);
  return res;
}

sexp sexp_ratio_compare (sexp ctx, sexp a, sexp b) {
  sexp_gc_var2(a2, b2);
  sexp_gc_preserve2(ctx, a2, b2);
  a2 = sexp_mul(ctx, sexp_ratio_numerator(a), sexp_ratio_denominator(b));
  b2 = sexp_mul(ctx, sexp_ratio_numerator(b), sexp_ratio_denominator(a));
  a2 = sexp_compare(ctx, a2, b2);
  sexp_gc_release2(ctx);
  return a2;
}

sexp sexp_ratio_round (sexp ctx, sexp a) {
  sexp_gc_var2(q, r);
  sexp_gc_preserve2(ctx, q, r);
  q = sexp_quotient(ctx, sexp_ratio_numerator(a), sexp_ratio_denominator(a));
  if ((sexp_ratio_denominator(a) == SEXP_TWO) && sexp_oddp(q)) {
    q = sexp_add(ctx, q, (sexp_exact_positivep(q) ? SEXP_ONE : SEXP_NEG_ONE));
  } else {
    r = sexp_remainder(ctx, sexp_ratio_numerator(a), sexp_ratio_denominator(a));
    r = sexp_mul(ctx, r, SEXP_TWO);
    if (sexp_exact_negativep(r)) {sexp_negate(r);}
    if (sexp_unbox_fixnum(sexp_compare(ctx, r, sexp_ratio_denominator(a))) > 0)
      q = sexp_add(ctx, q, (sexp_exact_negativep(sexp_ratio_numerator(a)) ? SEXP_NEG_ONE : SEXP_ONE));
  }
  sexp_gc_release2(ctx);
  return q;
}

sexp sexp_ratio_trunc (sexp ctx, sexp a) {
  return sexp_quotient(ctx, sexp_ratio_numerator(a), sexp_ratio_denominator(a));
}

sexp sexp_ratio_floor (sexp ctx, sexp a) {
  sexp_gc_var1(q);
  sexp_gc_preserve1(ctx, q);
  q = sexp_quotient(ctx, sexp_ratio_numerator(a), sexp_ratio_denominator(a));
  if (sexp_exact_negativep(sexp_ratio_numerator(a)))
    q = sexp_add(ctx, q, SEXP_NEG_ONE);
  sexp_gc_release1(ctx);
  return q;
}

sexp sexp_ratio_ceiling (sexp ctx, sexp a) {
  sexp_gc_var1(q);
  sexp_gc_preserve1(ctx, q);
  q = sexp_quotient(ctx, sexp_ratio_numerator(a), sexp_ratio_denominator(a));
  if (sexp_exact_positivep(sexp_ratio_numerator(a)))
    q = sexp_add(ctx, q, SEXP_ONE);
  sexp_gc_release1(ctx);
  return q;
}

#endif

double sexp_to_double (sexp ctx, sexp x) {
  if (sexp_flonump(x))
    return sexp_flonum_value(x);
  else if (sexp_fixnump(x))
    return sexp_fixnum_to_double(x);
  else if (sexp_bignump(x))
    return sexp_bignum_to_double(x);
#if SEXP_USE_RATIOS
  else if (sexp_ratiop(x))
    return sexp_ratio_to_double(ctx, x);
#endif
  else
    return 0.0;
}

/************************ complex numbers ****************************/

#if SEXP_USE_COMPLEX

static sexp sexp_complex_copy (sexp ctx, sexp a) {
  sexp_gc_var1(res);
  sexp_gc_preserve1(ctx, res);
  res = sexp_make_complex(ctx, sexp_complex_real(a), sexp_complex_imag(a));
  if (sexp_flonump(sexp_complex_real(a)))
    sexp_complex_real(a) = sexp_make_flonum(ctx, sexp_flonum_value(sexp_complex_real(a)));
  else if (sexp_bignump(sexp_complex_real(a)))
    sexp_complex_real(a) = sexp_copy_bignum(ctx, NULL, sexp_complex_real(a), 0);
  if (sexp_flonump(sexp_complex_imag(a)))
    sexp_complex_imag(a) = sexp_make_flonum(ctx, sexp_flonum_value(sexp_complex_imag(a)));
  else if (sexp_bignump(sexp_complex_imag(a)))
    sexp_complex_imag(a) = sexp_copy_bignum(ctx, NULL, sexp_complex_imag(a), 0);
  sexp_gc_release1(ctx);
  return res;
}

sexp sexp_complex_add (sexp ctx, sexp a, sexp b) {
  sexp_gc_var3(res, real, imag);
  sexp_gc_preserve3(ctx, res, real, imag);
  real = sexp_add(ctx, sexp_complex_real(a), sexp_complex_real(b));
  imag = sexp_add(ctx, sexp_complex_imag(a), sexp_complex_imag(b));
  res = sexp_make_complex(ctx, real, imag);
  sexp_gc_release3(ctx);
  return sexp_complex_normalize(res);
}

sexp sexp_complex_sub (sexp ctx, sexp a, sexp b) {
  sexp_gc_var2(res, tmp);
  sexp_gc_preserve2(ctx, res, tmp);
  tmp = sexp_complex_copy(ctx, b);
  sexp_negate_maybe_ratio(sexp_complex_real(tmp));
  sexp_negate_maybe_ratio(sexp_complex_imag(tmp));
  res = sexp_complex_add(ctx, a, tmp);
  sexp_gc_release2(ctx);
  return res;
}

sexp sexp_complex_mul (sexp ctx, sexp a, sexp b) {
  sexp_gc_var3(res, real, imag);
  sexp_gc_preserve3(ctx, res, real, imag);
  real = sexp_mul(ctx, sexp_complex_real(a), sexp_complex_real(b));
  res = sexp_mul(ctx, sexp_complex_imag(a), sexp_complex_imag(b));
  real = sexp_sub(ctx, real, res);
  imag = sexp_mul(ctx, sexp_complex_real(a), sexp_complex_imag(b));
  res = sexp_mul(ctx, sexp_complex_imag(a), sexp_complex_real(b));
  imag = sexp_add(ctx, imag, res);
  res = sexp_make_complex(ctx, real, imag);
  sexp_gc_release3(ctx);
  return sexp_complex_normalize(res);
}

/* (a + bi)    (ac + bd)      (bc - ad)    */
/* -------- = -----------  + ----------- i */
/* (c + di)   (c^2 + d^2)    (c^2 + d^2)   */

sexp sexp_complex_div (sexp ctx, sexp a, sexp b) {
  sexp_gc_var4(res, real, imag, denom);
  sexp_gc_preserve4(ctx, res, real, imag, denom);
  /* c^2 + d^2 */
  denom = sexp_mul(ctx, sexp_complex_real(b), sexp_complex_real(b));
  res = sexp_mul(ctx, sexp_complex_imag(b), sexp_complex_imag(b));
  denom = sexp_add(ctx, denom, res);
  /* ac + bd */
  real = sexp_mul(ctx, sexp_complex_real(a), sexp_complex_real(b));
  res = sexp_mul(ctx, sexp_complex_imag(a), sexp_complex_imag(b));
  real = sexp_add(ctx, real, res);
  real = sexp_div(ctx, real, denom);
  /* bc - ad */
  imag = sexp_mul(ctx, sexp_complex_imag(a), sexp_complex_real(b));
  res = sexp_mul(ctx, sexp_complex_real(a), sexp_complex_imag(b));
  imag = sexp_sub(ctx, imag, res);
  imag = sexp_div(ctx, imag, denom);
  res = sexp_make_complex(ctx, real, imag);
  sexp_gc_release4(ctx);
  return sexp_complex_normalize(res);
}

static sexp sexp_to_complex (sexp ctx, sexp x) {
#if SEXP_USE_RATIOS
  sexp_gc_var1(tmp);
#endif
  if (sexp_flonump(x) || sexp_fixnump(x) || sexp_bignump(x)) {
    return sexp_make_complex(ctx, x, SEXP_ZERO);
#if SEXP_USE_RATIOS
  } else if (sexp_ratiop(x)) {
    sexp_gc_preserve1(ctx, tmp);
    tmp = sexp_make_complex(ctx, SEXP_ZERO, SEXP_ZERO);
    sexp_complex_real(tmp) = sexp_make_flonum(ctx, sexp_to_double(ctx, x));
    sexp_gc_release1(ctx);
    return tmp;
#endif
  } else {
    return x;
  }
}

sexp sexp_complex_exp (sexp ctx, sexp z) {
  double e2x = exp(sexp_to_double(ctx, sexp_complex_real(z))),
    y = sexp_to_double(ctx, sexp_complex_imag(z));
  sexp_gc_var1(res);
  sexp_gc_preserve1(ctx, res);
  res = sexp_make_complex(ctx, SEXP_ZERO, SEXP_ZERO);
  sexp_complex_real(res) = sexp_make_flonum(ctx, e2x*cos(y));
  sexp_complex_imag(res) = sexp_make_flonum(ctx, e2x*sin(y));
  sexp_gc_release1(ctx);
  return res;
}

sexp sexp_complex_log (sexp ctx, sexp z) {
  double x = sexp_to_double(ctx, sexp_complex_real(z)),
    y = sexp_to_double(ctx, sexp_complex_imag(z));
  sexp_gc_var1(res);
  sexp_gc_preserve1(ctx, res);
  res = sexp_make_complex(ctx, SEXP_ZERO, SEXP_ZERO);
  sexp_complex_real(res) = sexp_make_flonum(ctx, log(sqrt(x*x + y*y)));
  sexp_complex_imag(res) = sexp_make_flonum(ctx, atan2(y, x));
  sexp_gc_release1(ctx);
  return res;
}

sexp sexp_complex_expt (sexp ctx, sexp a, sexp b) {
  sexp_gc_var1(res);
  sexp_gc_preserve1(ctx, res);
  res = sexp_to_complex(ctx, a);
  res = sexp_complex_log(ctx, res);
  res = sexp_mul(ctx, b, res);
  res = sexp_complex_exp(ctx, res);
  sexp_gc_release1(ctx);
  return res;
}

#if SEXP_USE_MATH

sexp sexp_complex_sqrt (sexp ctx, sexp z) {
  double x = sexp_to_double(ctx, sexp_complex_real(z)),
    y = sexp_to_double(ctx, sexp_complex_imag(z)), r;
  sexp_gc_var1(res);
  sexp_gc_preserve1(ctx, res);
  r = sqrt(x*x + y*y);
  res = sexp_make_complex(ctx, SEXP_ZERO, SEXP_ZERO);
  sexp_complex_real(res) = sexp_make_flonum(ctx, sqrt((x+r)/2));
  sexp_complex_imag(res) = sexp_make_flonum(ctx, ((y<-0.0)?-1:1)*sqrt((-x+r)/2));
  sexp_gc_release1(ctx);
  return res;
}

sexp sexp_complex_sin (sexp ctx, sexp z) {
  double x = sexp_to_double(ctx, sexp_complex_real(z)),
    y = sexp_to_double(ctx, sexp_complex_imag(z));
  sexp_gc_var1(res);
  sexp_gc_preserve1(ctx, res);
  res = sexp_make_complex(ctx, SEXP_ZERO, SEXP_ZERO);
  sexp_complex_real(res) = sexp_make_flonum(ctx, sin(x)*cosh(y));
  sexp_complex_imag(res) = sexp_make_flonum(ctx, cos(x)*sinh(y));
  sexp_gc_release1(ctx);
  return res;
}

sexp sexp_complex_cos (sexp ctx, sexp z) {
  double x = sexp_to_double(ctx, sexp_complex_real(z)),
    y = sexp_to_double(ctx, sexp_complex_imag(z));
  sexp_gc_var1(res);
  sexp_gc_preserve1(ctx, res);
  res = sexp_make_complex(ctx, SEXP_ZERO, SEXP_ZERO);
  sexp_complex_real(res) = sexp_make_flonum(ctx, cos(x)*cosh(y));
  sexp_complex_imag(res) = sexp_make_flonum(ctx, -sin(x)*sinh(y));
  sexp_gc_release1(ctx);
  return res;
}

sexp sexp_complex_tan (sexp ctx, sexp z) {
  sexp res;
  sexp_gc_var2(sin, cos);
  sexp_gc_preserve2(ctx, sin, cos);
  sin = sexp_complex_sin(ctx, z);
  cos = sexp_complex_cos(ctx, z);
  res = sexp_complex_div(ctx, sin, cos);
  sexp_gc_release2(ctx);
  return res;
}

sexp sexp_complex_asin (sexp ctx, sexp z) {
  sexp_gc_var3(res, tmp, tmp2);
  sexp_gc_preserve3(ctx, res, tmp, tmp2);
  res = sexp_complex_mul(ctx, z, z);
  res = sexp_sub(ctx, SEXP_ONE, res);
  res = sexp_sqrt(ctx, NULL, 1, res);
  tmp = sexp_make_complex(ctx, SEXP_ZERO, SEXP_ZERO);
  sexp_complex_real(tmp) = sexp_mul(ctx, SEXP_NEG_ONE, sexp_complex_imag(z));
  sexp_complex_imag(tmp) = sexp_complex_real(z);
  res = sexp_add(ctx, tmp, res);
  res = sexp_log(ctx, NULL, 1, res);
  tmp = sexp_make_complex(ctx, SEXP_ZERO, SEXP_NEG_ONE);
  res = sexp_mul(ctx, res, tmp);
  sexp_gc_release3(ctx);
  return res;
}

sexp sexp_complex_acos (sexp ctx, sexp z) {
  sexp_gc_var2(res, tmp);
  sexp_gc_preserve2(ctx, res, tmp);
  res = sexp_complex_asin(ctx, z);
  tmp = sexp_make_complex(ctx, SEXP_ZERO, SEXP_ZERO);
  sexp_complex_real(tmp) = sexp_make_flonum(ctx, acos(-1)/2);
  res = sexp_sub(ctx, tmp, res);
  sexp_gc_release2(ctx);
  return res;
}

sexp sexp_complex_atan (sexp ctx, sexp z) {
  sexp_gc_var3(res, tmp1, tmp2);
  sexp_gc_preserve3(ctx, res, tmp1, tmp2);
  tmp1 = sexp_make_complex(ctx, SEXP_ZERO, SEXP_ONE);
  tmp1 = sexp_complex_mul(ctx, z, tmp1);
  res = sexp_make_complex(ctx, SEXP_ONE, SEXP_ZERO);
  res = sexp_complex_sub(ctx, res, tmp1);
  res = sexp_complex_log(ctx, res);
  tmp2 = sexp_make_complex(ctx, SEXP_ONE, SEXP_ZERO);
  tmp2 = sexp_complex_add(ctx, tmp2, tmp1);
  tmp2 = sexp_complex_log(ctx, tmp2);
  res = sexp_complex_sub(ctx, res, tmp2);
  tmp1 = sexp_make_complex(ctx, SEXP_ZERO, SEXP_ONE);
  sexp_complex_imag(tmp1) = sexp_make_flonum(ctx, 0.5);
  res = sexp_complex_mul(ctx, res, tmp1);
  sexp_gc_release3(ctx);
  return res;
}

#endif
#endif

/****************** generic arithmetic ************************/

enum sexp_number_types {
  SEXP_NUM_NOT = 0,
  SEXP_NUM_FIX,
  SEXP_NUM_FLO,
  SEXP_NUM_BIG,
#if SEXP_USE_RATIOS
  SEXP_NUM_RAT,
#endif
#if SEXP_USE_COMPLEX
  SEXP_NUM_CPX,
#endif
};

enum sexp_number_combs {
  SEXP_NUM_NOT_NOT = 0,
  SEXP_NUM_NOT_FIX,
  SEXP_NUM_NOT_FLO,
  SEXP_NUM_NOT_BIG,
#if SEXP_USE_RATIOS
  SEXP_NUM_NOT_RAT,
#endif
#if SEXP_USE_COMPLEX
  SEXP_NUM_NOT_CPX,
#endif
  SEXP_NUM_FIX_NOT,
  SEXP_NUM_FIX_FIX,
  SEXP_NUM_FIX_FLO,
  SEXP_NUM_FIX_BIG,
#if SEXP_USE_RATIOS
  SEXP_NUM_FIX_RAT,
#endif
#if SEXP_USE_COMPLEX
  SEXP_NUM_FIX_CPX,
#endif
  SEXP_NUM_FLO_NOT,
  SEXP_NUM_FLO_FIX,
  SEXP_NUM_FLO_FLO,
  SEXP_NUM_FLO_BIG,
#if SEXP_USE_RATIOS
  SEXP_NUM_FLO_RAT,
#endif
#if SEXP_USE_COMPLEX
  SEXP_NUM_FLO_CPX,
#endif
  SEXP_NUM_BIG_NOT,
  SEXP_NUM_BIG_FIX,
  SEXP_NUM_BIG_FLO,
  SEXP_NUM_BIG_BIG,
#if SEXP_USE_RATIOS
  SEXP_NUM_BIG_RAT,
#endif
#if SEXP_USE_COMPLEX
  SEXP_NUM_BIG_CPX,
#endif
#if SEXP_USE_RATIOS
  SEXP_NUM_RAT_NOT,
  SEXP_NUM_RAT_FIX,
  SEXP_NUM_RAT_FLO,
  SEXP_NUM_RAT_BIG,
  SEXP_NUM_RAT_RAT,
#if SEXP_USE_COMPLEX
  SEXP_NUM_RAT_CPX,
#endif
#endif
#if SEXP_USE_COMPLEX
  SEXP_NUM_CPX_NOT,
  SEXP_NUM_CPX_FIX,
  SEXP_NUM_CPX_FLO,
  SEXP_NUM_CPX_BIG,
#if SEXP_USE_RATIOS
  SEXP_NUM_CPX_RAT,
#endif
  SEXP_NUM_CPX_CPX,
#endif
};

static int sexp_number_types[] =
#if SEXP_USE_RATIOS && SEXP_USE_COMPLEX
  {0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 2, 3, 4, 5, 0, 0};
#else
#if SEXP_USE_RATIOS || SEXP_USE_COMPLEX
  {0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 2, 3, 4, 0, 0, 0};
#else
  {0, 0, 1, 0, 0, 0, 0, 0, 0, 0, 0, 2, 3, 0, 0, 0, 0};
#endif
#endif

#define SEXP_NUM_NUMBER_TYPES (4 + SEXP_USE_RATIOS + SEXP_USE_COMPLEX)

static int sexp_number_type (sexp a) {
  return sexp_pointerp(a) ?
    (sexp_pointer_tag(a)<(sizeof(sexp_number_types)/sizeof(sexp_number_types[0]))
     ? sexp_number_types[sexp_pointer_tag(a)] : 0)
#if SEXP_USE_IMMEDIATE_FLONUMS
    : sexp_flonump(a) ? 2
#endif
    : sexp_fixnump(a);
}

sexp sexp_add (sexp ctx, sexp a, sexp b) {
  sexp_sint_t sum;
  int at=sexp_number_type(a), bt=sexp_number_type(b), t;
  sexp r=SEXP_VOID;
  sexp_gc_var1(tmp);
  sexp_gc_preserve1(ctx, tmp);
  if (at > bt) {r=a; a=b; b=r; t=at; at=bt; bt=t;}
  switch ((at * SEXP_NUM_NUMBER_TYPES) + bt) {
  case SEXP_NUM_NOT_NOT: case SEXP_NUM_NOT_FIX:
  case SEXP_NUM_NOT_FLO: case SEXP_NUM_NOT_BIG:
#if SEXP_USE_RATIOS
  case SEXP_NUM_NOT_RAT:
#endif
#if SEXP_USE_COMPLEX
  case SEXP_NUM_NOT_CPX:
#endif
    r = sexp_type_exception(ctx, NULL, SEXP_NUMBER, a);
    break;
  case SEXP_NUM_FIX_FIX:
    sum = sexp_unbox_fixnum(a) + sexp_unbox_fixnum(b);
    if ((sum < SEXP_MIN_FIXNUM) || (sum > SEXP_MAX_FIXNUM))
      r = sexp_add(ctx, tmp=sexp_fixnum_to_bignum(ctx, a), b);
    else
      r = sexp_make_fixnum(sum);
    break;
  case SEXP_NUM_FIX_FLO:
    r = a == SEXP_ZERO ? b : sexp_make_flonum(ctx, sexp_fixnum_to_double(a)+sexp_flonum_value(b));
    break;
  case SEXP_NUM_FIX_BIG:
    r = sexp_bignum_normalize(sexp_bignum_add_fixnum(ctx, b, a));
    break;
  case SEXP_NUM_FLO_FLO:
    r = sexp_fp_add(ctx, a, b);
    break;
  case SEXP_NUM_FLO_BIG:
    r = sexp_make_flonum(ctx, sexp_flonum_value(a) + sexp_bignum_to_double(b));
    break;
  case SEXP_NUM_BIG_BIG:
    r = sexp_bignum_normalize(sexp_bignum_add(ctx, NULL, b, a));
    break;
#if SEXP_USE_RATIOS
  case SEXP_NUM_FLO_RAT:
    r = sexp_make_flonum(ctx, sexp_flonum_value(a) + sexp_ratio_to_double(ctx, b));
    break;
  case SEXP_NUM_FIX_RAT:
  case SEXP_NUM_BIG_RAT:
    a = tmp = sexp_make_ratio(ctx, a, SEXP_ONE);
    /* ... FALLTHROUGH ... */
  case SEXP_NUM_RAT_RAT:
    r = sexp_ratio_add(ctx, a, b);
    break;
#endif
#if SEXP_USE_COMPLEX
#if SEXP_USE_RATIOS
  case SEXP_NUM_RAT_CPX:
#endif
  case SEXP_NUM_FLO_CPX:
  case SEXP_NUM_FIX_CPX:
  case SEXP_NUM_BIG_CPX:
    a = tmp = sexp_make_complex(ctx, a, SEXP_ZERO);
    /* ... FALLTHROUGH ... */
  case SEXP_NUM_CPX_CPX:
    r = sexp_complex_add(ctx, a, b);
    break;
#endif
  }
  sexp_gc_release1(ctx);
  return r;
}

sexp sexp_sub (sexp ctx, sexp a, sexp b) {
#if SEXP_USE_FLONUMS
  int negatep=0;
#endif
  int at=sexp_number_type(a), bt=sexp_number_type(b);
  sexp r=SEXP_VOID;
  sexp_gc_var2(tmp1, tmp2);
  sexp_gc_preserve2(ctx, tmp1, tmp2);
  switch ((at * SEXP_NUM_NUMBER_TYPES) + bt) {
  case SEXP_NUM_NOT_NOT: case SEXP_NUM_NOT_FIX:
  case SEXP_NUM_NOT_FLO: case SEXP_NUM_NOT_BIG:
#if SEXP_USE_RATIOS
  case SEXP_NUM_NOT_RAT:
#endif
#if SEXP_USE_COMPLEX
  case SEXP_NUM_NOT_CPX:
#endif
    r = sexp_type_exception(ctx, NULL, SEXP_NUMBER, a);
    break;
  case SEXP_NUM_FIX_NOT: case SEXP_NUM_FLO_NOT: case SEXP_NUM_BIG_NOT:
#if SEXP_USE_RATIOS
  case SEXP_NUM_RAT_NOT:
#endif
#if SEXP_USE_COMPLEX
  case SEXP_NUM_CPX_NOT:
#endif
    r = sexp_type_exception(ctx, NULL, SEXP_NUMBER, b);
    break;
  case SEXP_NUM_FIX_FIX:
    r = sexp_fx_sub(a, b);      /* VM catches this case */
    break;
  case SEXP_NUM_FIX_FLO:
    r = sexp_make_flonum(ctx, a==SEXP_ZERO ? -sexp_flonum_value(b) : sexp_fixnum_to_double(a)-sexp_flonum_value(b));
    break;
  case SEXP_NUM_FIX_BIG:
    tmp1 = sexp_fixnum_to_bignum(ctx, a);
    r = sexp_bignum_sub(ctx, NULL, b, tmp1);
    sexp_negate_exact(r);
    r = sexp_bignum_normalize(r);
    break;
  case SEXP_NUM_FLO_FIX:
    r = sexp_make_flonum(ctx, sexp_flonum_value(a) - sexp_fixnum_to_double(b));
    break;
  case SEXP_NUM_FLO_FLO:
    r = sexp_fp_sub(ctx, a, b);
    break;
  case SEXP_NUM_FLO_BIG:
    r = sexp_make_flonum(ctx, sexp_flonum_value(a) - sexp_bignum_to_double(b));
    break;
  case SEXP_NUM_BIG_FIX:
    tmp1 = sexp_fixnum_to_bignum(ctx, b);
    r = sexp_bignum_normalize(sexp_bignum_sub(ctx, NULL, a, tmp1));
    break;
  case SEXP_NUM_BIG_FLO:
    r = sexp_make_flonum(ctx, sexp_bignum_to_double(a) - sexp_flonum_value(b));
    break;
  case SEXP_NUM_BIG_BIG:
    r = sexp_bignum_normalize(sexp_bignum_sub(ctx, NULL, a, b));
    break;
#if SEXP_USE_RATIOS
  case SEXP_NUM_FLO_RAT:
    r = sexp_make_flonum(ctx, sexp_flonum_value(a) - sexp_ratio_to_double(ctx, b));
    break;
  case SEXP_NUM_RAT_FLO:
    r = sexp_make_flonum(ctx, sexp_ratio_to_double(ctx, a) - sexp_flonum_value(b));
    break;
  case SEXP_NUM_RAT_FIX:
  case SEXP_NUM_RAT_BIG:
    tmp1 = a; a = b; b = tmp1;
    negatep = 1;
    /* ... FALLTHROUGH ... */
  case SEXP_NUM_FIX_RAT:
  case SEXP_NUM_BIG_RAT:
    a = tmp1 = sexp_make_ratio(ctx, a, SEXP_ONE);
    /* ... FALLTHROUGH ... */
  case SEXP_NUM_RAT_RAT:
    tmp2 = sexp_make_ratio(ctx, sexp_ratio_numerator(b), sexp_ratio_denominator(b));
    if (sexp_bignump(sexp_ratio_numerator(tmp2)))
      sexp_ratio_numerator(tmp2) = sexp_copy_bignum(ctx, NULL, sexp_ratio_numerator(tmp2), 0);
    sexp_negate_exact(sexp_ratio_numerator(tmp2));
    r = sexp_ratio_add(ctx, a, tmp2);
    if (negatep) {
      sexp_negate_maybe_ratio(r);
    }
    break;
#endif
#if SEXP_USE_COMPLEX
#if SEXP_USE_RATIOS
  case SEXP_NUM_RAT_CPX:
    a = tmp1 = sexp_make_flonum(ctx, sexp_ratio_to_double(ctx, a));
    goto complex_sub;
  case SEXP_NUM_CPX_RAT:
    b = tmp1 = sexp_make_flonum(ctx, sexp_ratio_to_double(ctx, b));
    /* ... FALLTHROUGH ... */
#endif
  case SEXP_NUM_CPX_FLO:
  case SEXP_NUM_CPX_FIX:
  case SEXP_NUM_CPX_BIG:
    tmp1 = a; a = b; b = tmp1;
    negatep = 1;
    /* ... FALLTHROUGH ... */
  case SEXP_NUM_FLO_CPX:
  case SEXP_NUM_FIX_CPX:
  case SEXP_NUM_BIG_CPX:
    a = tmp1 = sexp_make_complex(ctx, a, SEXP_ZERO);
    /* ... FALLTHROUGH ... */
  case SEXP_NUM_CPX_CPX:
#if SEXP_USE_RATIOS
  complex_sub:
#endif
    r = sexp_complex_sub(ctx, a, b);
    if (negatep) {
      if (sexp_complexp(r)) {
        r = sexp_complex_copy(ctx, r);
        sexp_negate_maybe_ratio(sexp_complex_real(r));
        sexp_negate_maybe_ratio(sexp_complex_imag(r));
      } else {
        sexp_negate_maybe_ratio(r);
      }
    }
    break;
#endif
  }
  sexp_gc_release2(ctx);
  return r;
}

sexp sexp_mul (sexp ctx, sexp a, sexp b) {
  sexp_lsint_t prod;
  int at=sexp_number_type(a), bt=sexp_number_type(b), t;
  sexp r=SEXP_VOID;
  sexp_gc_var1(tmp);
  sexp_gc_preserve1(ctx, tmp);
  if (at > bt) {r=a; a=b; b=r; t=at; at=bt; bt=t;}
  switch ((at * SEXP_NUM_NUMBER_TYPES) + bt) {
  case SEXP_NUM_NOT_NOT: case SEXP_NUM_NOT_FIX:
  case SEXP_NUM_NOT_FLO: case SEXP_NUM_NOT_BIG:
#if SEXP_USE_RATIOS
  case SEXP_NUM_NOT_RAT:
#endif
    r = sexp_type_exception(ctx, NULL, SEXP_NUMBER, a);
    break;
  case SEXP_NUM_FIX_FIX:
    prod = lsint_mul_sint(lsint_from_sint(sexp_unbox_fixnum(a)), sexp_unbox_fixnum(b));
    if (!lsint_is_fixnum(prod))
      r = sexp_mul(ctx, tmp=sexp_fixnum_to_bignum(ctx, a), b);
    else
      r = sexp_make_fixnum(lsint_to_sint(prod));
    break;
  case SEXP_NUM_FIX_FLO:
    r = (a==SEXP_ZERO ? a : sexp_make_flonum(ctx, sexp_fixnum_to_double(a)*sexp_flonum_value(b)));
    break;
  case SEXP_NUM_FIX_BIG:
    r = sexp_bignum_fxmul(ctx, NULL, b, sexp_unbox_fx_abs(a), 0);
    sexp_bignum_sign(r) = sexp_fx_sign(a) * sexp_bignum_sign(b);
    r = sexp_bignum_normalize(r);
    break;
  case SEXP_NUM_FLO_FLO:
    r = sexp_fp_mul(ctx, a, b);
    break;
  case SEXP_NUM_FLO_BIG:
    r = sexp_make_flonum(ctx, sexp_flonum_value(a) * sexp_bignum_to_double(b));
    break;
  case SEXP_NUM_BIG_BIG:
    r = sexp_bignum_normalize(sexp_bignum_mul(ctx, NULL, a, b));
    break;
#if SEXP_USE_RATIOS
  case SEXP_NUM_FLO_RAT:
    r = sexp_make_flonum(ctx, sexp_flonum_value(a) * sexp_ratio_to_double(ctx, b));
    break;
  case SEXP_NUM_FIX_RAT:
  case SEXP_NUM_BIG_RAT:
    a = tmp = sexp_make_ratio(ctx, a, SEXP_ONE);
    /* ... FALLTHROUGH ... */
  case SEXP_NUM_RAT_RAT:
    r = sexp_ratio_mul(ctx, a, b);
    break;
#endif
#if SEXP_USE_COMPLEX
#if SEXP_USE_RATIOS
  case SEXP_NUM_RAT_CPX:
#endif
  case SEXP_NUM_FLO_CPX:
  case SEXP_NUM_FIX_CPX:
  case SEXP_NUM_BIG_CPX:
    a = tmp = sexp_make_complex(ctx, a, SEXP_ZERO);
    /* ... FALLTHROUGH ... */
  case SEXP_NUM_CPX_CPX:
    r = sexp_complex_mul(ctx, a, b);
    break;
#endif
  }
  sexp_gc_release1(ctx);
  return r;
}

sexp sexp_div (sexp ctx, sexp a, sexp b) {
  int at=sexp_number_type(a), bt=sexp_number_type(b);
#if ! SEXP_USE_RATIOS
  double f;
#endif
  sexp r=SEXP_VOID;
  sexp_gc_var1(tmp);
  sexp_gc_preserve1(ctx, tmp);
  switch ((at * SEXP_NUM_NUMBER_TYPES) + bt) {
  case SEXP_NUM_NOT_NOT: case SEXP_NUM_NOT_FIX:
  case SEXP_NUM_NOT_FLO: case SEXP_NUM_NOT_BIG:
#if SEXP_USE_RATIOS
  case SEXP_NUM_NOT_RAT:
#endif
    r = sexp_type_exception(ctx, NULL, SEXP_NUMBER, a);
    break;
  case SEXP_NUM_FIX_NOT: case SEXP_NUM_FLO_NOT: case SEXP_NUM_BIG_NOT:
#if SEXP_USE_RATIOS
  case SEXP_NUM_RAT_NOT:
#endif
    r = sexp_type_exception(ctx, NULL, SEXP_NUMBER, b);
    break;
  case SEXP_NUM_FIX_FIX:
#if SEXP_USE_RATIOS
    tmp = sexp_make_ratio(ctx, a, b);
    r = sexp_ratio_normalize(ctx, tmp, SEXP_FALSE);
#else
    f = sexp_fixnum_to_double(a) / sexp_fixnum_to_double(b);
    r = ((f == trunc(f)) ? sexp_make_fixnum((sexp_sint_t)f)
         : sexp_make_flonum(ctx, f));
#endif
    break;
  case SEXP_NUM_FIX_FLO:
    r = sexp_make_flonum(ctx, sexp_fixnum_to_double(a)/sexp_flonum_value(b));
    break;
  case SEXP_NUM_FIX_BIG:
#if SEXP_USE_RATIOS
    tmp = sexp_make_ratio(ctx, a, b);
    r = sexp_ratio_normalize(ctx, tmp, SEXP_FALSE);
#else
    r = sexp_make_flonum(ctx, sexp_fixnum_to_double(a)/sexp_bignum_to_double(b));
#endif
    break;
  case SEXP_NUM_FLO_FIX:
    r = sexp_make_flonum(ctx, sexp_flonum_value(a)/sexp_fixnum_to_double(b));
    break;
  case SEXP_NUM_FLO_FLO:
    r = sexp_fp_div(ctx, a, b);
    break;
  case SEXP_NUM_FLO_BIG:
    r = sexp_make_flonum(ctx, sexp_flonum_value(a) / sexp_bignum_to_double(b));
    break;
  case SEXP_NUM_BIG_FIX:
#if SEXP_USE_RATIOS
    tmp = sexp_make_ratio(ctx, a, b);
    r = sexp_ratio_normalize(ctx, tmp, SEXP_FALSE);
    break;
#else
    b = tmp = sexp_fixnum_to_bignum(ctx, b);
#endif
    /* ... FALLTHROUGH if ! SEXP_USE_RATIOS ... */
  case SEXP_NUM_BIG_BIG:
#if SEXP_USE_RATIOS
    tmp = sexp_make_ratio(ctx, a, b);
    r = sexp_ratio_normalize(ctx, tmp, SEXP_FALSE);
#else
    r = sexp_bignum_quot_rem(ctx, &tmp, a, b);
    if (sexp_bignum_normalize(tmp) != SEXP_ZERO)
      r = sexp_make_flonum(ctx, sexp_bignum_to_double(a)
                           / sexp_bignum_to_double(b));
    else
      r = sexp_bignum_normalize(r);
#endif
    break;
  case SEXP_NUM_BIG_FLO:
    r = sexp_make_flonum(ctx, sexp_bignum_to_double(a) / sexp_flonum_value(b));
    break;
#if SEXP_USE_RATIOS
  case SEXP_NUM_FLO_RAT:
    r = sexp_make_flonum(ctx, sexp_flonum_value(a) / sexp_ratio_to_double(ctx, b));
    break;
  case SEXP_NUM_RAT_FLO:
    r = sexp_make_flonum(ctx, sexp_ratio_to_double(ctx, a) / sexp_flonum_value(b));
    break;
  case SEXP_NUM_RAT_FIX:
  case SEXP_NUM_RAT_BIG:
    b = tmp = sexp_make_ratio(ctx, b, SEXP_ONE);
    /* ... FALLTHROUGH ... */
  case SEXP_NUM_FIX_RAT:
  case SEXP_NUM_BIG_RAT:
    if (!sexp_ratiop(a))
      a = tmp = sexp_make_ratio(ctx, a, SEXP_ONE);
    /* ... FALLTHROUGH ... */
  case SEXP_NUM_RAT_RAT:
    r = sexp_ratio_div(ctx, a, b);
    break;
#endif
#if SEXP_USE_COMPLEX
#if SEXP_USE_RATIOS
  case SEXP_NUM_CPX_RAT:
    b = tmp = sexp_make_flonum(ctx, sexp_ratio_to_double(ctx, b));
    /* ... FALLTHROUGH ... */
#endif
  case SEXP_NUM_CPX_FLO:
  case SEXP_NUM_CPX_FIX:
  case SEXP_NUM_CPX_BIG:
    b = tmp = sexp_make_complex(ctx, b, SEXP_ZERO);
    /* ... FALLTHROUGH ... */
#if SEXP_USE_RATIOS
  case SEXP_NUM_RAT_CPX:
    if (sexp_ratiop(a))
      a = tmp = sexp_make_flonum(ctx, sexp_ratio_to_double(ctx, a));
    /* ... FALLTHROUGH ... */
#endif
  case SEXP_NUM_FLO_CPX:
  case SEXP_NUM_FIX_CPX:
  case SEXP_NUM_BIG_CPX:
    if (!sexp_complexp(a))
      a = tmp = sexp_make_complex(ctx, a, SEXP_ZERO);
    /* ... FALLTHROUGH ... */
  case SEXP_NUM_CPX_CPX:
    r = sexp_complex_div(ctx, a, b);
    break;
#endif
  }
  sexp_gc_release1(ctx);
  return r;
}

sexp sexp_to_inexact (sexp ctx, sexp a) {
  if (sexp_fixnump(a)) return sexp_fixnum_to_flonum(ctx, a);
  if (sexp_bignump(a)) return sexp_make_flonum(ctx, sexp_bignum_to_double(a));
  return a;
}

sexp sexp_quotient (sexp ctx, sexp a, sexp b) {
  int at=sexp_number_type(a), bt=sexp_number_type(b);
  sexp r=SEXP_VOID;
  sexp_gc_var1(tmp);
  if (b == SEXP_ONE) return a;
  sexp_gc_preserve1(ctx, tmp);
  switch ((at * SEXP_NUM_NUMBER_TYPES) + bt) {
  case SEXP_NUM_NOT_NOT: case SEXP_NUM_NOT_FIX:
  case SEXP_NUM_NOT_FLO: case SEXP_NUM_NOT_BIG:
    r = sexp_type_exception(ctx, NULL, SEXP_FIXNUM, a);
    break;
  case SEXP_NUM_FIX_NOT: case SEXP_NUM_FLO_NOT: case SEXP_NUM_BIG_NOT:
    r = sexp_type_exception(ctx, NULL, SEXP_FIXNUM, b);
    break;
  case SEXP_NUM_FLO_FIX: case SEXP_NUM_FLO_FLO: case SEXP_NUM_FLO_BIG:
#if SEXP_USE_RATIOS
  case SEXP_NUM_FLO_RAT:
#endif
    if (sexp_flonum_value(a) != trunc(sexp_flonum_value(a))) {
      r = sexp_type_exception(ctx, NULL, SEXP_FIXNUM, a);
    } else {
      tmp = sexp_bignum_normalize(sexp_double_to_bignum(ctx, sexp_flonum_value(a)));
      tmp = sexp_quotient(ctx, tmp, b);
      r = sexp_to_inexact(ctx, tmp);
    }
    break;
#if SEXP_USE_RATIOS
  case SEXP_NUM_RAT_FIX: case SEXP_NUM_RAT_BIG: case SEXP_NUM_RAT_RAT:
#endif
#if SEXP_USE_COMPLEX
  case SEXP_NUM_FLO_CPX: case SEXP_NUM_CPX_FIX: case SEXP_NUM_CPX_FLO:
  case SEXP_NUM_CPX_BIG: case SEXP_NUM_CPX_CPX:
#if SEXP_USE_RATIOS
  case SEXP_NUM_CPX_RAT:
#endif
#endif
    r = sexp_type_exception(ctx, NULL, SEXP_FIXNUM, a);
    break;
  case SEXP_NUM_FIX_FLO: case SEXP_NUM_BIG_FLO:
#if SEXP_USE_RATIOS
  case SEXP_NUM_RAT_FLO:
#endif
    if (sexp_flonum_value(b) != trunc(sexp_flonum_value(b))) {
      r = sexp_type_exception(ctx, NULL, SEXP_FIXNUM, b);
    } else {
      tmp = sexp_bignum_normalize(sexp_double_to_bignum(ctx, sexp_flonum_value(b)));
      tmp = sexp_quotient(ctx, a, tmp);
      r = sexp_to_inexact(ctx, tmp);
    }
    break;
#if SEXP_USE_RATIOS
  case SEXP_NUM_FIX_RAT: case SEXP_NUM_BIG_RAT:
#endif
#if SEXP_USE_COMPLEX
  case SEXP_NUM_FIX_CPX: case SEXP_NUM_BIG_CPX:
#endif
    r = sexp_type_exception(ctx, NULL, SEXP_FIXNUM, b);
    break;
  case SEXP_NUM_FIX_FIX:
    r = sexp_fx_div(a, b);
    if ((sexp_sint_t)a < 0 && (sexp_sint_t)b < 0 && (sexp_sint_t)r < 0) {
      r = sexp_quotient(ctx, tmp=sexp_fixnum_to_bignum(ctx, a), b);
    }
    break;
  case SEXP_NUM_FIX_BIG:
    r = SEXP_ZERO;
    break;
  case SEXP_NUM_BIG_FIX:
    b = tmp = sexp_fixnum_to_bignum(ctx, b);
    /* ... FALLTHROUGH ... */
  case SEXP_NUM_BIG_BIG:
    r = sexp_bignum_normalize(sexp_bignum_quotient(ctx, a, b));
    break;
  }
  sexp_gc_release1(ctx);
  return r;
}

sexp sexp_remainder (sexp ctx, sexp a, sexp b) {
  int at=sexp_number_type(a), bt=sexp_number_type(b);
  sexp r=SEXP_VOID;
  sexp_gc_var1(tmp);
  if (b == SEXP_ONE) return SEXP_ZERO;
  sexp_gc_preserve1(ctx, tmp);
  switch ((at * SEXP_NUM_NUMBER_TYPES) + bt) {
  case SEXP_NUM_NOT_NOT: case SEXP_NUM_NOT_FIX:
  case SEXP_NUM_NOT_FLO: case SEXP_NUM_NOT_BIG:
    r = sexp_type_exception(ctx, NULL, SEXP_FIXNUM, a);
    break;
  case SEXP_NUM_FIX_NOT: case SEXP_NUM_FLO_NOT: case SEXP_NUM_BIG_NOT:
    r = sexp_type_exception(ctx, NULL, SEXP_FIXNUM, b);
    break;
  case SEXP_NUM_FLO_FIX: case SEXP_NUM_FLO_FLO: case SEXP_NUM_FLO_BIG:
#if SEXP_USE_RATIOS
  case SEXP_NUM_FLO_RAT:
#endif
    if (isinf(sexp_flonum_value(a)) ||
        sexp_flonum_value(a) != trunc(sexp_flonum_value(a))) {
      r = sexp_type_exception(ctx, NULL, SEXP_FIXNUM, a);
    } else if (bt == SEXP_NUM_FLO && isinf(sexp_flonum_value(b))) {
      r = sexp_type_exception(ctx, NULL, SEXP_FIXNUM, b);
    } else {
      tmp = sexp_bignum_normalize(sexp_double_to_bignum(ctx, sexp_flonum_value(a)));
      tmp = sexp_remainder(ctx, tmp, b);
      r = sexp_to_inexact(ctx, tmp);
    }
    break;
#if SEXP_USE_RATIOS
  case SEXP_NUM_RAT_FIX: case SEXP_NUM_RAT_BIG: case SEXP_NUM_RAT_RAT:
#endif
#if SEXP_USE_COMPLEX
  case SEXP_NUM_FLO_CPX: case SEXP_NUM_CPX_FIX: case SEXP_NUM_CPX_FLO:
  case SEXP_NUM_CPX_BIG: case SEXP_NUM_CPX_CPX:
#if SEXP_USE_RATIOS
  case SEXP_NUM_CPX_RAT:
#endif
#endif
    r = sexp_type_exception(ctx, NULL, SEXP_FIXNUM, a);
    break;
  case SEXP_NUM_FIX_FLO: case SEXP_NUM_BIG_FLO:
#if SEXP_USE_RATIOS
  case SEXP_NUM_RAT_FLO:
#endif
    if (isinf(sexp_flonum_value(b)) ||
        sexp_flonum_value(b) != trunc(sexp_flonum_value(b))) {
      r = sexp_type_exception(ctx, NULL, SEXP_FIXNUM, b);
    } else {
      tmp = sexp_bignum_normalize(sexp_double_to_bignum(ctx, sexp_flonum_value(b)));
      tmp = sexp_remainder(ctx, a, tmp);
      r = sexp_to_inexact(ctx, tmp);
    }
    break;
#if SEXP_USE_RATIOS
  case SEXP_NUM_FIX_RAT: case SEXP_NUM_BIG_RAT:
#endif
#if SEXP_USE_COMPLEX
  case SEXP_NUM_FIX_CPX: case SEXP_NUM_BIG_CPX:
#endif
    r = sexp_type_exception(ctx, NULL, SEXP_FIXNUM, b);
    break;
  case SEXP_NUM_FIX_FIX:
    r = sexp_fx_rem(a, b);
    break;
  case SEXP_NUM_FIX_BIG:
    r = a;
    break;
  case SEXP_NUM_BIG_FIX:
    r = sexp_bignum_fxrem(ctx, a, sexp_unbox_fixnum(b));
    break;
  case SEXP_NUM_BIG_BIG:
    r = sexp_bignum_normalize(sexp_bignum_remainder(ctx, a, b));
    break;
  }
  sexp_gc_release1(ctx);
  return r;
}

sexp sexp_compare (sexp ctx, sexp a, sexp b) {
  int at=sexp_number_type(a), bt=sexp_number_type(b);
  sexp r=SEXP_VOID;
  double f, g;
  sexp_gc_var1(tmp);
  sexp_gc_preserve1(ctx, tmp);
  if (at > bt) {
    r = sexp_compare(ctx, b, a);
    if (!sexp_exceptionp(r)) { sexp_negate(r); }
  } else {
    switch ((at * SEXP_NUM_NUMBER_TYPES) + bt) {
    case SEXP_NUM_NOT_NOT: case SEXP_NUM_NOT_FIX:
    case SEXP_NUM_NOT_FLO: case SEXP_NUM_NOT_BIG:
#if SEXP_USE_COMPLEX
    case SEXP_NUM_CPX_CPX: case SEXP_NUM_FIX_CPX:
    case SEXP_NUM_FLO_CPX: case SEXP_NUM_BIG_CPX:
#if SEXP_USE_RATIOS
    case SEXP_NUM_RAT_CPX:
#endif
#endif
      r = sexp_type_exception(ctx, NULL, SEXP_NUMBER, a);
      break;
    case SEXP_NUM_FIX_FIX:
      r = sexp_make_fixnum(sexp_unbox_fixnum(a) - sexp_unbox_fixnum(b));
      break;
    case SEXP_NUM_FIX_FLO:
      if (isinf(sexp_flonum_value(b))) {
        r = sexp_flonum_value(b) > 0 ? SEXP_NEG_ONE : SEXP_ONE;
      } else if (isnan(sexp_flonum_value(b))) {
        r = sexp_xtype_exception(ctx, NULL, "can't compare NaN", b);
      } else {
        r = sexp_compare(ctx, a, tmp=sexp_inexact_to_exact(ctx, NULL, 1, b));
      }
      break;
    case SEXP_NUM_FIX_BIG:
      if ((sexp_bignum_hi(b) > 1) ||
          (sexp_bignum_data(b)[0] > SEXP_MAX_FIXNUM))
        r = sexp_make_fixnum(sexp_bignum_sign(b) < 0 ? 1 : -1);
      else
        r = sexp_make_fixnum(sexp_unbox_fixnum(a) - (sexp_sint_t)sexp_bignum_data(b)[0]);
      break;
    case SEXP_NUM_FLO_FLO:
      f = sexp_flonum_value(a);
      g = sexp_flonum_value(b);
      if (isnan(f))
        r = sexp_xtype_exception(ctx, NULL, "can't compare NaN", a);
      else if (isnan(g))
        r = sexp_xtype_exception(ctx, NULL, "can't compare NaN", b);
      else
        r = sexp_make_fixnum(f < g ? -1 : f == g ? 0 : 1);
      break;
    case SEXP_NUM_FLO_BIG:
      f = sexp_flonum_value(a);
      if (isinf(f)) {
        r = f > 0 ? SEXP_ONE : SEXP_NEG_ONE;
        break;
      } else if (isnan(f)) {
        r = sexp_xtype_exception(ctx, NULL, "can't compare NaN", a);
        break;
      } else {
        a = tmp = sexp_double_to_bignum(ctx, f);
      }
      /* ... FALLTHROUGH ... */
    case SEXP_NUM_BIG_BIG:
      r = sexp_make_fixnum(sexp_bignum_compare(a, b));
      break;
#if SEXP_USE_RATIOS
    case SEXP_NUM_FLO_RAT:
      f = sexp_flonum_value(a);
      if (isinf(f)) {
        r = f > 0 ? SEXP_ONE : SEXP_NEG_ONE;
      } else if (isnan(f)) {
        r = sexp_xtype_exception(ctx, NULL, "can't compare NaN", a);
      } else {
        r = sexp_compare(ctx, tmp=sexp_inexact_to_exact(ctx, NULL, 1, a), b);
      }
      break;
    case SEXP_NUM_FIX_RAT:
    case SEXP_NUM_BIG_RAT:
      a = tmp = sexp_make_ratio(ctx, a, SEXP_ONE);
      /* ... FALLTHROUGH ... */
    case SEXP_NUM_RAT_RAT:
      r = sexp_ratio_compare(ctx, a, b);
      break;
#endif
    default:
      r = sexp_xtype_exception(ctx, NULL, "unknown comparison", a);
      break;
    }
  }
  sexp_gc_release1(ctx);
  return r;
}

#endif