More portable bignums that don't have to rely on gcc 128bit arithmetics extension

SEXP_USE_CUSTOM_LONG_LONGS currently needs SEXP_64_BIT

AUTHORS

@@ -34,6 +34,7 @@ Thanks to the following people for patches and bug reports:
   * Bakul Shah
   * Ben Mather
   * Ben Weaver
+  * Bertrand Augereau
   * Bruno Deferrari
   * Doug Currie
   * Derrick Eddington

bignum.c

@@ -35,35 +35,44 @@ sexp sexp_fixnum_to_bignum (sexp ctx, sexp a) {
   return res;
 }
 
-sexp sexp_make_integer (sexp ctx, sexp_lsint_t x) {
+sexp sexp_make_integer_from_lsint (sexp ctx, sexp_lsint_t x) {
   sexp res;
-  if ((SEXP_MIN_FIXNUM <= x) && (x <= SEXP_MAX_FIXNUM)) {
-    res = sexp_make_fixnum(x);
+  if (lsint_is_fixnum(x)) {
+    res = sexp_make_fixnum(lsint_to_sint(x));
   } else {
     res = sexp_make_bignum(ctx, 1);
-    if (x < 0) {
+    if (lsint_lt_0(x)) {
       sexp_bignum_sign(res) = -1;
-      sexp_bignum_data(res)[0] = (sexp_uint_t)-x;
+      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)x;
+      sexp_bignum_data(res)[0] = (sexp_uint_t)lsint_to_sint(x);
     }
   }
   return res;
 }
 
-sexp sexp_make_unsigned_integer (sexp ctx, sexp_luint_t x) {
+sexp sexp_make_unsigned_integer_from_luint (sexp ctx, sexp_luint_t x) {
   sexp res;
-  if (x <= SEXP_MAX_FIXNUM) {
-    res = sexp_make_fixnum(x);
+  if (luint_is_fixnum(x)) {
+    res = sexp_make_fixnum(luint_to_uint(x));
   } else {
     res = sexp_make_bignum(ctx, 1);
     sexp_bignum_sign(res) = 1;
-    sexp_bignum_data(res)[0] = (sexp_uint_t)x;
+    sexp_bignum_data(res)[0] = luint_to_uint(x);
   }
   return res;
 }
 
+#if !SEXP_USE_CUSTOM_LONG_LONGS
+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  /* !SEXP_USE_CUSTOM_LONG_LONGS */
+
 #define double_trunc_10s_digit(f) (trunc((f)/10.0)*10.0)
 #define double_10s_digit(f) ((f)-double_trunc_10s_digit(f))
 
@@ -200,9 +209,9 @@ sexp sexp_bignum_fxmul (sexp ctx, sexp d, sexp a, sexp_uint_t b, int offset) {
   tmp = d;
   data = sexp_bignum_data(d);
   for (i=0; i<len; i++) {
-    n = (sexp_luint_t)adata[i]*b + carry;
-    data[i+offset] = (sexp_uint_t)n;
-    carry = n >> (sizeof(sexp_uint_t)*8);
+    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)
@@ -216,13 +225,13 @@ sexp sexp_bignum_fxmul (sexp ctx, sexp d, sexp a, sexp_uint_t b, int offset) {
 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 = 0;
+  sexp_luint_t n = luint_from_uint(0);
   for (i=len-1; i>=offset; i--) {
-    n = (n << sizeof(sexp_uint_t)*8) + data[i];
-    q = (sexp_uint_t)(n / b);
-    r = (sexp_uint_t)(n - (sexp_luint_t)q * b);
+    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 = r;
+    n = luint_from_uint(r);
   }
   return r;
 }
@@ -230,7 +239,7 @@ sexp_uint_t sexp_bignum_fxdiv (sexp ctx, sexp a, sexp_uint_t b, int offset) {
 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 = 0;
+  sexp_luint_t n = luint_from_uint(0);
   if (b > 0) {
     q = b - 1;
     if ((b & q) == 0)
@@ -241,11 +250,11 @@ sexp sexp_bignum_fxrem (sexp ctx, sexp a, sexp_sint_t b) {
     return sexp_xtype_exception(ctx, NULL, "divide by zero", a);
   }
   for (i=len-1; i>=0; i--) {
-    n = (n << sizeof(sexp_uint_t)*8) + data[i];
-    q = (sexp_uint_t)(n / b0);
-    n -= (sexp_luint_t)q * b0;
+    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)n);
+  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,
@@ -526,38 +535,38 @@ sexp sexp_bignum_quot_rem (sexp ctx, sexp *rem, sexp a, sexp b) {
     sexp_bignum_data(x)[off] = 0;
     if (off > 0) sexp_bignum_data(x)[off-1] = 0;
     off = alen - blen + 1;
-    dn = (((sexp_luint_t)sexp_bignum_data(a1)[alen-1]
-           << (sizeof(sexp_uint_t)*8))
-          + sexp_bignum_data(a1)[alen-2]);
-    dd = (((sexp_luint_t)sexp_bignum_data(b1)[blen-1]
-           << (sizeof(sexp_uint_t)*8))
-          + sexp_bignum_data(b1)[blen-2]);
+    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 &&
-        sexp_bignum_data(a1)[alen-1] < ((sexp_luint_t)1<<(sizeof(sexp_uint_t)*4)) &&
-        sexp_bignum_data(b1)[blen-1] < ((sexp_luint_t)1<<(sizeof(sexp_uint_t)*4))) {
-      dn = (dn << (sizeof(sexp_uint_t)*4))
-        + (sexp_bignum_data(a1)[alen-3] >> (sizeof(sexp_uint_t)*4));
-      dd = (dd << (sizeof(sexp_uint_t)*4))
-        + (sexp_bignum_data(b1)[blen-3] >> (sizeof(sexp_uint_t)*4));
+        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 = dn / dd;
-    if (d == 0) {
-      dn = (((sexp_luint_t)sexp_bignum_data(a1)[alen-1]
-             << (sizeof(sexp_uint_t)*8))
-            + sexp_bignum_data(a1)[alen-2]);
-      dd = sexp_bignum_data(b1)[blen-1];
-      if (sexp_bignum_data(a1)[alen-1] < ((sexp_luint_t)1<<(sizeof(sexp_uint_t)*4)) &&
-          sexp_bignum_data(b1)[blen-1] < ((sexp_luint_t)1<<(sizeof(sexp_uint_t)*4))) {
-        dn = (dn << (sizeof(sexp_uint_t)*4))
-          + (sexp_bignum_data(a1)[alen-3] >> (sizeof(sexp_uint_t)*4));
-        dd = (dd << (sizeof(sexp_uint_t)*4))
-          + (sexp_bignum_data(b1)[blen-2] >> (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 = dn / dd;
+      d = luint_div(dn, dd);
       off--;
     }
-    dhi = d >> (sizeof(sexp_uint_t)*8);
-    dlo = d & (((sexp_luint_t)1<<(sizeof(sexp_uint_t)*8))-1);
+    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 */
@@ -1430,11 +1439,11 @@ sexp sexp_mul (sexp ctx, sexp a, sexp b) {
     r = sexp_type_exception(ctx, NULL, SEXP_NUMBER, a);
     break;
   case SEXP_NUM_FIX_FIX:
-    prod = (sexp_lsint_t)sexp_unbox_fixnum(a) * sexp_unbox_fixnum(b);
-    if ((prod < SEXP_MIN_FIXNUM) || (prod > SEXP_MAX_FIXNUM))
+    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(prod);
+      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)));

include/chibi/bignum.h

@@ -1,5 +1,5 @@
 /*  bignum.h -- header for bignum utilities                   */
-/*  Copyright (c) 2009-2015 Alex Shinn.  All rights reserved. */
+/*  Copyright (c) 2009-2018 Alex Shinn.  All rights reserved. */
 /*  BSD-style license: http://synthcode.com/license.txt       */
 
 #ifndef SEXP_BIGNUM_H
@@ -7,7 +7,19 @@
 
 #include "chibi/eval.h"
 
-#if (SEXP_64_BIT) && defined(__GNUC__)
+#if SEXP_USE_CUSTOM_LONG_LONGS
+#include <stdint.h>
+typedef struct
+{
+  uint64_t hi;
+  uint64_t lo;
+} sexp_luint_t;
+typedef struct
+{
+  int64_t  hi;
+  uint64_t lo;
+} sexp_lsint_t;
+#elif SEXP_64_BIT
 typedef unsigned int uint128_t __attribute__((mode(TI)));
 typedef int sint128_t __attribute__((mode(TI)));
 typedef uint128_t sexp_luint_t;
@@ -17,6 +29,337 @@ typedef unsigned long long sexp_luint_t;
 typedef long long sexp_lsint_t;
 #endif
 
+#if !SEXP_USE_CUSTOM_LONG_LONGS
+
+#define lsint_from_sint(v) ((sexp_lsint_t)v)
+#define luint_from_uint(v) ((sexp_luint_t)v)
+#define lsint_to_sint(v) ((sexp_sint_t)v)
+#define luint_to_uint(v) ((sexp_uint_t)v)
+#define lsint_negate(v) (-((sexp_lsint_t)v))
+#define luint_eq(a, b) (((sexp_luint_t)a)==((sexp_luint_t)b))
+#define luint_lt(a, b) (((sexp_luint_t)a)<((sexp_luint_t)b))
+#define lsint_lt_0(a) (((sexp_lsint_t)a)<0)
+#define luint_shl(a, shift) (((sexp_luint_t)a)<<(shift))
+#define luint_shr(a, shift) (((sexp_luint_t)a)>>(shift))
+#define luint_add(a, b) (((sexp_luint_t)a)+((sexp_luint_t)b))
+#define luint_add_uint(a, b) (((sexp_luint_t)a)+((sexp_uint_t)b))
+#define luint_sub(a, b) (((sexp_luint_t)a)-((sexp_luint_t)b))
+#define luint_mul_uint(a, b) (((sexp_luint_t)a)*((sexp_uint_t)b))
+#define lsint_mul_sint(a, b) (((sexp_lsint_t)a)*((sexp_sint_t)b))
+#define luint_div(a, b) (((sexp_luint_t)a)/((sexp_luint_t)b))
+#define luint_div_uint(a, b) (((sexp_luint_t)a)/((sexp_luint_t)b))
+#define luint_and(a, b) (((sexp_luint_t)a)&((sexp_luint_t)b))
+#define luint_is_fixnum(x) (((sexp_luint_t)x)<=SEXP_MAX_FIXNUM)
+#define lsint_is_fixnum(x) ((SEXP_MIN_FIXNUM <= ((sexp_lsint_t)x)) && (((sexp_lsint_t)x) <= SEXP_MAX_FIXNUM))
+
+#else
+
+static inline sexp_luint_t luint_from_lsint(sexp_lsint_t v) {
+  sexp_luint_t result;
+  result.hi = v.hi;
+  result.lo = v.lo;
+  return result;
+}
+
+static inline sexp_lsint_t lsint_from_luint(sexp_luint_t v) {
+  sexp_lsint_t result;
+  result.hi = v.hi;
+  result.lo = v.lo;
+  return result;
+}
+
+static inline sexp_lsint_t lsint_from_sint(sexp_sint_t v) {
+  sexp_lsint_t result;
+  result.hi = v >> 63;
+  result.lo = v;
+  return result;
+}
+
+static inline sexp_luint_t luint_from_uint(sexp_uint_t v) {
+  sexp_luint_t result;
+  result.hi = 0;
+  result.lo = v;
+  return result;
+}
+
+static inline sexp_sint_t lsint_to_sint(sexp_lsint_t v) {
+  return v.lo;
+}
+
+static inline sexp_uint_t luint_to_uint(sexp_luint_t v) {
+  return v.lo;
+}
+
+static inline sexp_lsint_t lsint_negate(sexp_lsint_t v) {
+  sexp_luint_t a;
+  a.hi = ~v.hi;
+  a.lo = ~v.lo;
+
+  uint64_t aLoLo = a.lo & 0xFFFFFFFF;
+  uint64_t aLoHi = a.lo >> 32;
+  uint64_t aHiLo = a.hi & 0xFFFFFFFF;
+  uint64_t aHiHi = a.hi >> 32;
+ 
+  uint64_t carry;
+  uint64_t sumLoLo = aLoLo + 1;
+  carry = sumLoLo >> 32;
+  uint64_t resultLoLo = sumLoLo & 0xFFFFFFFF;
+
+  uint64_t sumLoHi = aLoHi + carry;
+  uint64_t resultLoHi = sumLoHi & 0xFFFFFFFF;
+  carry = sumLoHi >> 32;
+
+  uint64_t sumHiLo = aHiLo + carry;
+  uint64_t resultHiLo = sumHiLo & 0xFFFFFFFF;
+  carry = sumHiLo >> 32;
+
+  uint64_t sumHiHi = aHiHi + carry;
+  uint64_t resultHiHi = sumHiHi & 0xFFFFFFFF;
+  //carry = sumHiHi >> 32;
+
+  sexp_lsint_t result;
+  result.hi = (resultHiHi << 32) | resultHiLo;
+  result.lo = (resultLoHi << 32) | resultLoLo;
+  return result;
+}
+
+static inline int luint_eq(sexp_luint_t a, sexp_luint_t b) {
+  return (a.hi == b.hi) && (a.lo == b.lo);
+}
+
+static inline int luint_lt(sexp_luint_t a, sexp_luint_t b) {
+  if (a.hi < b.hi)
+    return 1;
+  else if (a.hi > b.hi)
+    return 0;
+  else
+    return a.lo < b.lo;
+}
+
+static inline int lsint_lt_0(sexp_lsint_t a) {
+  return a.hi < 0;
+}
+
+static inline sexp_luint_t luint_shl(sexp_luint_t v, size_t shift) {
+  if (shift == 0)
+    return v;
+  sexp_luint_t result;
+  if (shift >= 64) {
+    result.hi = v.lo << (shift - 64);
+    result.lo = 0;
+  } else {
+    result.hi = (v.hi << shift) | (v.lo >> (64-shift));
+    result.lo = v.lo << shift;
+  }
+  return result;
+}
+
+static inline sexp_luint_t luint_shr(sexp_luint_t v, size_t shift) {
+  if (shift == 0)
+    return v;
+  sexp_luint_t result;
+  if (shift >= 64) {
+    result.hi = 0;
+    result.lo = v.hi >> (shift - 64);
+  } else {
+    result.hi = v.hi >> shift;
+    result.lo = (v.lo >> shift) | (v.hi << (64-shift));
+  }
+  return result;
+}
+
+static inline sexp_luint_t luint_add(sexp_luint_t a, sexp_luint_t b) {
+  uint64_t aLoLo = a.lo & 0xFFFFFFFF;
+  uint64_t aLoHi = a.lo >> 32;
+  uint64_t aHiLo = a.hi & 0xFFFFFFFF;
+  uint64_t aHiHi = a.hi >> 32;
+  uint64_t bLoLo = b.lo & 0xFFFFFFFF;
+  uint64_t bLoHi = b.lo >> 32;
+  uint64_t bHiLo = b.hi & 0xFFFFFFFF;
+  uint64_t bHiHi = b.hi >> 32;
+
+  uint64_t carry;
+  uint64_t sumLoLo = (aLoLo + bLoLo);
+  carry = sumLoLo >> 32;
+  uint64_t resultLoLo = sumLoLo & 0xFFFFFFFF;
+
+  uint64_t sumLoHi = (aLoHi + bLoHi) + carry;
+  uint64_t resultLoHi = sumLoHi & 0xFFFFFFFF;
+  carry = sumLoHi >> 32;
+
+  uint64_t sumHiLo = (aHiLo + bHiLo) + carry;
+  uint64_t resultHiLo = sumHiLo & 0xFFFFFFFF;
+  carry = sumHiLo >> 32;
+
+  uint64_t sumHiHi = (aHiHi + bHiHi) + carry;
+  uint64_t resultHiHi = sumHiHi & 0xFFFFFFFF;
+  //carry = sumHiHi >> 32;
+
+  sexp_luint_t result;
+  result.hi = (resultHiHi << 32) | resultHiLo;
+  result.lo = (resultLoHi << 32) | resultLoLo;
+  return result;
+}
+
+static inline sexp_luint_t luint_add_uint(sexp_luint_t a, sexp_uint_t b) {
+  uint64_t aLoLo = a.lo & 0xFFFFFFFF;
+  uint64_t aLoHi = a.lo >> 32;
+  uint64_t aHiLo = a.hi & 0xFFFFFFFF;
+  uint64_t aHiHi = a.hi >> 32;
+  uint64_t bLoLo = b & 0xFFFFFFFF;
+  uint64_t bLoHi = b >> 32;
+
+  uint64_t carry;
+  uint64_t sumLoLo = (aLoLo + bLoLo);
+  carry = sumLoLo >> 32;
+  uint64_t resultLoLo = sumLoLo & 0xFFFFFFFF;
+
+  uint64_t sumLoHi = (aLoHi + bLoHi) + carry;
+  uint64_t resultLoHi = sumLoHi & 0xFFFFFFFF;
+  carry = sumLoHi >> 32;
+
+  uint64_t sumHiLo = aHiLo + carry;
+  uint64_t resultHiLo = sumHiLo & 0xFFFFFFFF;
+  carry = sumHiLo >> 32;
+
+  uint64_t sumHiHi = aHiHi + carry;
+  uint64_t resultHiHi = sumHiHi & 0xFFFFFFFF;
+  //carry = sumHiHi >> 32;
+
+  sexp_luint_t result;
+  result.hi = (resultHiHi << 32) | resultHiLo;
+  result.lo = (resultLoHi << 32) | resultLoLo;
+  return result;
+}
+
+static inline sexp_luint_t luint_sub(sexp_luint_t a, sexp_luint_t b) {
+  sexp_luint_t negB;
+  negB.hi = ~b.hi;
+  negB.lo = ~b.lo;
+  return luint_add(a, luint_add_uint(negB, 1));
+}
+
+static inline sexp_luint_t luint_mul_uint(sexp_luint_t a, sexp_uint_t b) {
+  uint64_t aLoLo = a.lo & 0xFFFFFFFF;
+  uint64_t aLoHi = a.lo >> 32;
+  uint64_t aHiLo = a.hi & 0xFFFFFFFF;
+  uint64_t aHiHi = a.hi >> 32;
+
+  uint64_t bLo = b & 0xFFFFFFFF;
+  uint64_t bHi = b >> 32;
+
+  sexp_luint_t resultBLo, resultBHi;
+  {
+    sexp_luint_t prodLoLo;
+    prodLoLo.hi = 0;
+    prodLoLo.lo = aLoLo * bLo;
+
+    sexp_luint_t prodLoHi;
+    prodLoHi.hi = (aLoHi * bLo) >> 32;
+    prodLoHi.lo = (aLoHi * bLo) << 32;
+
+    sexp_luint_t prodHiLo;
+    prodHiLo.hi = aHiLo * bLo;
+    prodHiLo.lo = 0;
+
+    sexp_luint_t prodHiHi;
+    prodHiHi.hi = (aHiHi * bLo) << 32;
+    prodHiHi.lo = 0;
+
+    resultBLo = luint_add(luint_add(luint_add(prodLoLo, prodLoHi), prodHiLo), prodHiHi);
+  }
+  {
+    sexp_luint_t prodLoLo;
+    prodLoLo.hi = 0;
+    prodLoLo.lo = aLoLo * bHi;
+
+    sexp_luint_t prodLoHi;
+    prodLoHi.hi = (aLoHi * bHi) >> 32;
+    prodLoHi.lo = (aLoHi * bHi) << 32;
+
+    sexp_luint_t prodHiLo;
+    prodHiLo.hi = aHiLo * bHi;
+    prodHiLo.lo = 0;
+
+    sexp_luint_t prodHiHi;
+    prodHiHi.hi = (aHiHi * bHi) << 32;
+    prodHiHi.lo = 0;
+
+    resultBHi = luint_add(luint_add(luint_add(prodLoLo, prodLoHi), prodHiLo), prodHiHi);
+  }
+
+  sexp_luint_t result = luint_add(resultBLo, luint_shl(resultBHi, 32));
+  
+  return result;
+}
+
+static inline sexp_lsint_t lsint_mul_sint(sexp_lsint_t a, sexp_sint_t b) {
+  if (lsint_lt_0(a)) {
+    sexp_luint_t minusA = luint_from_lsint(lsint_negate(a));
+    if (b < 0)
+      return lsint_from_luint(luint_mul_uint(minusA, (sexp_uint_t)-b));
+    else
+      return lsint_negate(lsint_from_luint(luint_mul_uint(minusA, (sexp_uint_t)b)));
+  } else {
+    if (b < 0)
+      return lsint_negate(lsint_from_luint(luint_mul_uint(luint_from_lsint(a), (sexp_uint_t)-b)));
+    else
+      return lsint_from_luint(luint_mul_uint(luint_from_lsint(a), (sexp_uint_t)b));
+  }
+}
+
+static inline sexp_luint_t luint_div(sexp_luint_t a, sexp_luint_t b) {
+  if (luint_lt(a, b))
+    return luint_from_uint(0);
+  else if (luint_eq(a, b))
+    return luint_from_uint(1);
+
+  sexp_luint_t quotient = luint_from_uint(0);
+  sexp_luint_t remainder = luint_from_uint(0);
+
+  for (int i = 0; i < 128; i++) {
+    quotient = luint_shl(quotient, 1);
+
+    remainder = luint_shl(remainder, 1);
+    remainder.lo |= (a.hi >> 63) & 1;
+    a = luint_shl(a, 1);
+
+    if (!(luint_lt(remainder, b))) {
+      remainder = luint_sub(remainder, b);
+      quotient.lo |= 1;
+    }
+  }
+
+  return quotient;
+}
+
+static inline sexp_luint_t luint_div_uint(sexp_luint_t a, sexp_uint_t b) {
+  return luint_div(a, luint_from_uint(b));
+}
+
+static inline sexp_luint_t luint_and(sexp_luint_t a, sexp_luint_t b) {
+  sexp_luint_t result;
+  result.hi = a.hi & b.hi;
+  result.lo = a.lo & b.lo;
+  return result;
+}
+
+static inline int luint_is_fixnum(sexp_luint_t x) {
+  return (x.hi == 0) && (x.lo <= SEXP_MAX_FIXNUM);
+}
+
+static inline int lsint_is_fixnum(sexp_lsint_t x) {
+  if (x.hi > 0)
+      return 0;
+  else if (x.hi == 0)
+      return x.lo <= SEXP_MAX_FIXNUM;
+  else if (x.hi == -1)
+      return SEXP_MIN_FIXNUM <= x.lo;
+  else return 0;
+}
+
+#endif
+
 SEXP_API sexp_sint_t sexp_bignum_compare (sexp a, sexp b);
 SEXP_API sexp sexp_compare (sexp ctx, sexp a, sexp b);
 SEXP_API sexp sexp_make_bignum (sexp ctx, sexp_uint_t len);

include/chibi/features.h

@@ -301,6 +301,12 @@
 #endif
 #endif
 
+/* for bignum support, need a double long to store long*long */
+/* gcc supports uint128_t, otherwise we need a custom struct */
+#ifndef SEXP_USE_CUSTOM_LONG_LONGS
+#define SEXP_USE_CUSTOM_LONG_LONGS (SEXP_64_BIT && !defined(__GNUC__))
+#endif
+
 #ifndef SEXP_USE_NO_FEATURES
 #define SEXP_USE_NO_FEATURES 0
 #endif

include/chibi/sexp.h

@@ -891,8 +891,15 @@ SEXP_API int sexp_idp(sexp x);
 #endif
 
 #if SEXP_USE_BIGNUMS
+SEXP_API sexp sexp_make_integer_from_lsint(sexp ctx, sexp_lsint_t x);
+SEXP_API sexp sexp_make_unsigned_integer_from_luint(sexp ctx, sexp_luint_t x);
+#if SEXP_USE_CUSTOM_LONG_LONGS
+#define sexp_make_integer(ctx, x) sexp_make_fixnum(x)
+#define sexp_make_unsigned_integer(ctx, x) sexp_make_fixnum(x)
+#else
 SEXP_API sexp sexp_make_integer(sexp ctx, sexp_lsint_t x);
 SEXP_API sexp sexp_make_unsigned_integer(sexp ctx, sexp_luint_t x);
+#endif
 #define sexp_exact_integerp(x) (sexp_fixnump(x) || sexp_bignump(x))
 #else
 #define sexp_make_integer(ctx, x) sexp_make_fixnum(x)
@@ -1532,7 +1539,7 @@ SEXP_API sexp sexp_read_symbol (sexp ctx, sexp in, int init, int internp);
 SEXP_API sexp sexp_read_number (sexp ctx, sexp in, int base, int exactp);
 #if SEXP_USE_BIGNUMS
 SEXP_API sexp sexp_read_bignum (sexp ctx, sexp in, sexp_uint_t init,
-				signed char sign, sexp_uint_t base);
+                                signed char sign, sexp_uint_t base);
 SEXP_API sexp sexp_write_bignum (sexp ctx, sexp a, sexp out, sexp_uint_t base);
 #endif
 SEXP_API sexp sexp_read_float_tail(sexp ctx, sexp in, double whole, int negp);

tests/build/build-opts.txt

@@ -15,6 +15,7 @@ CPPFLAGS=-DSEXP_USE_IMMEDIATE_FLONUMS=1
 CPPFLAGS=-DSEXP_USE_BIGNUMS=0
 CPPFLAGS=-DSEXP_USE_COMPLEX=0
 CPPFLAGS=-DSEXP_USE_RATIOS=0
+CPPFLAGS=-DSEXP_USE_CUSTOM_LONG_LONGS=1
 CPPFLAGS=-DSEXP_USE_MATH=0
 CPPFLAGS=-DSEXP_WARN_UNDEFS=0
 CPPFLAGS=-DSEXP_USE_HUFF_SYMS=0

vm.c

@@ -1729,11 +1729,11 @@ sexp sexp_apply (sexp ctx, sexp proc, sexp args) {
     sexp_context_top(ctx) = --top;
 #if SEXP_USE_BIGNUMS
     if (sexp_fixnump(tmp1) && sexp_fixnump(tmp2)) {
-      prod = (sexp_lsint_t)sexp_unbox_fixnum(tmp1) * sexp_unbox_fixnum(tmp2);
-      if ((prod < SEXP_MIN_FIXNUM) || (prod > SEXP_MAX_FIXNUM))
+      prod = lsint_mul_sint(lsint_from_sint(sexp_unbox_fixnum(tmp1)), sexp_unbox_fixnum(tmp2));
+      if (!lsint_is_fixnum(prod))
         _ARG1 = sexp_mul(ctx, tmp1=sexp_fixnum_to_bignum(ctx, tmp1), tmp2);
       else
-        _ARG1 = sexp_make_fixnum(prod);
+        _ARG1 = sexp_make_fixnum(lsint_to_sint(prod));
     }
     else {
       _ARG1 = sexp_mul(ctx, tmp1, tmp2);