fix precision in exact->inexact on bignums

The strategy of converting a bignum to a flonum by converting on word
boundaries can lose one bit of precision. (If the use of a word
boundary causes a single bit to get rounded away, but the first bit of
the next word is non-zero, then the rounding might have been down when
it should have been up.)

Avoid the problem by aligning relative to the high bit, instead.

pkgs/racket-test-core/tests/racket/number.rktl

@@ -3196,6 +3196,96 @@
   
   (test #t list? (filter n-digit-has-nth-root? (build-list 5000 (lambda (x) (+ x 1))))))
 
+;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+;; exact->inexact precision on bignums (round-trip and proper rounding)
+
+(define max-53-bit-number (sub1 (arithmetic-shift 1 53)))
+
+(define (check-conversion 53-bit-number)
+  ;; Any 53-bit integer fits in a 64-bit floating point:
+  (unless (= 53-bit-number (inexact->exact (exact->inexact 53-bit-number)))
+    (error 'random-exact->inexact "round-trip failed ~s" 53-bit-number))
+  
+  ;; The same holds if we shift by up to (- 1023 52):
+  (define (check-shift p)
+    (define n2 (arithmetic-shift 53-bit-number p))
+    (unless (= n2 (inexact->exact (exact->inexact n2)))
+      (error 'random-exact->inexact "round-trip of shifted failed ~s" n2)))
+  (check-shift (- 1023 52))
+  (for ([i 10])
+    (check-shift (random (- 1023 52))))
+  
+  ;; The same holds if we shift by up to (- -1022 52):
+  (define (check-div p)
+    (define n2 (/ 53-bit-number (arithmetic-shift 1 (- p))))
+    (unless (= n2 (inexact->exact (exact->inexact n2)))
+      (error 'random-exact->inexact "round-trip of shifted failed ~s" n2)))
+  (check-div (- (+ 1022 52)))
+  (for ([i 10])
+    (check-div (- (random (+ 1022 52)))))
+  
+  ;; Helper for checking rounding:
+  (define (check-random-pairs check-shift-pair)
+    (check-shift-pair 1 0)
+    (check-shift-pair (- 1023 52 1) 0)
+    (check-shift-pair 1 (- 1023 52 2))
+    (for ([i 10])
+      (define zeros (add1 (random (- 1023 52 3))))
+      (define extra (random (- 1023 52 1 zeros)))
+      (check-shift-pair zeros extra)))
+  
+  ;; If we add a zero bit and then a non-zero bit anywhere later,
+  ;; conversion to inexact should round down.
+  (define (check-shift-plus-bits-to-truncate num-zeros extra-p)
+    (define n2 (arithmetic-shift
+                (bitwise-ior (arithmetic-shift 53-bit-number (add1 num-zeros))
+                             1)
+                extra-p))
+    (define n3 (inexact->exact (exact->inexact n2)))
+    (unless (= n3 (arithmetic-shift 53-bit-number (+ num-zeros 1 extra-p)))
+      (error 'random-exact->inexact "truncating round failed ~s" n2)))
+  (check-random-pairs check-shift-plus-bits-to-truncate)
+  
+  ;; If we add a one bit and then a non-zero bit anywhere later,
+  ;; conversion to inexact should round up.
+  (unless (= 53-bit-number max-53-bit-number)
+    (define (check-shift-plus-bits-to-up num-one-then-zeros extra-p)
+      (define n2 (arithmetic-shift
+                  (bitwise-ior (arithmetic-shift 
+                                (bitwise-ior (arithmetic-shift 53-bit-number 1)
+                                             1)
+                                num-one-then-zeros)
+                               1)
+                  extra-p))
+      (define n3 (inexact->exact (exact->inexact n2)))
+      (unless (= n3 (arithmetic-shift (add1 53-bit-number) (+ num-one-then-zeros 1 extra-p)))
+        (error 'random-exact->inexact "round up failed ~s" n2)))
+    (check-random-pairs check-shift-plus-bits-to-up))
+
+  ;; If we add a one bit and then only zero bits,
+  ;; conversion to inexact should round to even.
+  (unless (= 53-bit-number max-53-bit-number)
+    (define (check-shift-plus-bits-to-even num-one-then-zeros extra-p)
+      (define n2 (arithmetic-shift 
+                  (bitwise-ior (arithmetic-shift 53-bit-number 1)
+                               1)
+                  (+ num-one-then-zeros extra-p)))
+      (define n3 (inexact->exact (exact->inexact n2)))
+      (unless (= n3 (arithmetic-shift (if (even? 53-bit-number)
+                                          53-bit-number
+                                          (add1 53-bit-number))
+                                      (+ num-one-then-zeros 1 extra-p)))
+        (error 'random-exact->inexact "round to even failed ~s" n2)))
+    (check-random-pairs check-shift-plus-bits-to-even)))
+  
+(check-conversion max-53-bit-number)
+(for ([i 100])
+  (check-conversion
+   ;; Random 53-bit number:
+   (+ (arithmetic-shift 1 52)
+      (arithmetic-shift (random (expt 2 24)) 24)
+      (random (expt 2 28)))))
+
 ;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 ;; exact->inexact precision (thanks to Neil Toronto)
 

racket/collects/compiler/private/xform.rkt

@@ -903,7 +903,7 @@
 	       _isnan __isfinited __isnanl __isnan
                __isinff __isinfl isnanf isinff __isinfd __isnanf __isnand __isinf
                __inline_isnanl __inline_isnan
-               __builtin_popcount
+               __builtin_popcount __builtin_clz
                _Generic
                __inline_isinff __inline_isinfl __inline_isinfd __inline_isnanf __inline_isnand __inline_isinf
                floor floorl ceil ceill round roundl fmod fmodl modf modfl fabs fabsl __maskrune _errno __errno

racket/src/configure

@@ -5793,6 +5793,32 @@ if test "${has_builtin_popcount}" = "yes" ; then
 
 $as_echo "#define MZ_HAS_BUILTIN_POPCOUNT 1" >>confdefs.h
 
+fi
+
+ msg="for __builtin_clz"
+{ $as_echo "$as_me:${as_lineno-$LINENO}: checking $msg" >&5
+$as_echo_n "checking $msg... " >&6; }
+cat confdefs.h - <<_ACEOF >conftest.$ac_ext
+/* end confdefs.h.  */
+
+    int main(int argc, char **argv) {
+      unsigned int i = argc;
+      return __builtin_clz(i);
+    }
+_ACEOF
+if ac_fn_c_try_link "$LINENO"; then :
+  has_builtin_clz=yes
+else
+  has_builtin_clz=no
+fi
+rm -f core conftest.err conftest.$ac_objext \
+    conftest$ac_exeext conftest.$ac_ext
+{ $as_echo "$as_me:${as_lineno-$LINENO}: result: $has_builtin_clz" >&5
+$as_echo "$has_builtin_clz" >&6; }
+if test "${has_builtin_clz}" = "yes" ; then
+
+$as_echo "#define MZ_HAS_BUILTIN_CLZ 1" >>confdefs.h
+
 fi
 
 if test "${enable_backtrace}" = "yes" ; then

racket/src/racket/configure.ac

@@ -1276,6 +1276,18 @@ if test "${has_builtin_popcount}" = "yes" ; then
    AC_DEFINE(MZ_HAS_BUILTIN_POPCOUNT,1,[Has __builtin_popcount])
 fi
 
+[ msg="for __builtin_clz" ]
+AC_MSG_CHECKING($msg)
+AC_LINK_IFELSE([AC_LANG_SOURCE([
+    int main(int argc, char **argv) {
+      unsigned int i = argc;
+      return __builtin_clz(i);
+    }])], has_builtin_clz=yes, has_builtin_clz=no)
+AC_MSG_RESULT($has_builtin_clz)
+if test "${has_builtin_clz}" = "yes" ; then
+   AC_DEFINE(MZ_HAS_BUILTIN_CLZ,1,[Has __builtin_clz])
+fi
+
 if test "${enable_backtrace}" = "yes" ; then
   GC2OPTIONS="$GC2OPTIONS -DMZ_GC_BACKTRACE"
 fi

racket/src/racket/mzconfig.h.in

@@ -74,6 +74,9 @@ typedef unsigned long uintptr_t;
 /* When __builtin_popcount() is available: */
 #undef MZ_HAS_BUILTIN_POPCOUNT
 
+/* When __builtin_clz() is available: */
+#undef MZ_HAS_BUILTIN_CLZ
+
 /* Enable futures: */
 #undef MZ_USE_FUTURES
 

racket/src/racket/src/bgnfloat.inc

@@ -27,7 +27,7 @@ FP_TYPE SCHEME_BIGNUM_TO_FLOAT_INFO(const Scheme_Object *n, intptr_t skip, intpt
   FP_TYPE d;
 
   nl = SCHEME_BIGLEN(n);
-  na = SCHEME_BIGDIG(n) + nl;
+  na = SCHEME_BIGDIG(n);
 
   skipped = nl;
 
@@ -38,14 +38,63 @@ FP_TYPE SCHEME_BIGNUM_TO_FLOAT_INFO(const Scheme_Object *n, intptr_t skip, intpt
       return FP_scheme_floating_point_nzero;
   } else
     nl -= skip;
-	
-  d = FP_ZEROx;
-  while (nl--) {
-    d = FP_TYPE_MULT(d, FP_TYPE_FROM_DOUBLE(BIG_RADIX));
-    d = FP_TYPE_PLUS(d, FP_TYPE_FROM_UINTPTR(*(--na)));
-    if (IS_FLOAT_INF(d))
-      break;
-    --skipped;
+
+  if (!nl)
+    d = FP_ZEROx;
+  else if (nl == 1) {
+    d = FP_TYPE_FROM_UINTPTR(*na);
+    skipped = 0;
+  } else {
+    /* We'll get all the bits that matter in the first word or two,
+       and we won't lose precision as long as we shift so that the
+       highest bit in a word is non-zero */
+    bigdig b = na[nl-1];
+    int delta;
+
+    delta = mz_clz(b);
+    if (delta) {
+      /* zero bits in the highest word => pull in bits from the
+         second-highest word */
+      b = (b << delta) + (na[nl-2] >> (WORD_SIZE - delta));
+    }
+    if (sizeof(FP_TYPE) <= sizeof(bigdig)) {
+      /* one bigdig is enough, and the last bit is certainly
+         not needed, but it needs to summarize whether there
+         are any more non-zero bits in the number */
+      if (!(b & 0x1) && any_nonzero_digits(na, nl-1, delta))
+        b |= 0x1;
+      d = FP_TYPE_FROM_UINTPTR(b);
+    } else {
+      /* Need to look at a second word, possibly pulling in bits from
+         a third word */
+      d = FP_TYPE_FROM_UINTPTR(b);
+      d = FP_TYPE_MULT(d, FP_TYPE_FROM_DOUBLE(BIG_RADIX));
+      b = (na[nl-2] << delta);
+      if ((nl > 2) && delta)
+        b += (na[nl-3] >> (WORD_SIZE - delta));
+      if (!(b & 0x1) && (nl > 2) && any_nonzero_digits(na, nl-2, delta))
+        b |= 0x1;
+      d = FP_TYPE_PLUS(d, FP_TYPE_FROM_UINTPTR(b));
+      d = FP_TYPE_DIV(d, FP_TYPE_FROM_DOUBLE(BIG_RADIX));
+    }
+    /* Shift `d` back down by delta: */
+    if (delta)
+      d = FP_TYPE_DIV(d, FP_TYPE_POW(FP_TYPE_FROM_DOUBLE(2.0),
+                                     FP_TYPE_FROM_INT(delta)));
+    nl--;
+
+    /* Shift `d` up by remaining bignum words */
+    if (_skipped) {
+      while (nl--) {
+        d = FP_TYPE_MULT(d, FP_TYPE_FROM_DOUBLE(BIG_RADIX));
+        if (IS_FLOAT_INF(d))
+          break;
+        --skipped;
+      }
+    } else {
+      d = FP_TYPE_MULT(d, FP_TYPE_POW(FP_TYPE_FROM_DOUBLE(2.0),
+                                      FP_TYPE_FROM_UINTPTR(nl * WORD_SIZE)));
+    }
   }
 
   if (_skipped)
@@ -151,6 +200,7 @@ Scheme_Object *SCHEME_BIGNUM_FROM_FLOAT(FP_TYPE d)
 #undef FP_TYPE_MULT
 #undef FP_TYPE_PLUS
 #undef FP_TYPE_DIV
+#undef FP_TYPE_POW
 #undef FP_TYPE_FROM_INT
 #undef FP_TYPE_GREATER_OR_EQV
 #undef FP_TYPE_MINUS

racket/src/racket/src/bignum.c

@@ -1422,6 +1422,48 @@ static void bignum_add1_inplace(Scheme_Object **_stk_o)
     *_stk_o = bignum_copy(*_stk_o, carry);
 }
 
+XFORM_NONGCING static int mz_clz(uintptr_t n)
+{
+#ifdef MZ_HAS_BUILTIN_CLZ
+# if defined(SIXTY_FOUR_BIT_INTEGERS) || defined(USE_LONG_LONG_FOR_BIGDIG)
+  uintptr_t hi = (n >> (WORD_SIZE >> 1));
+  if (hi)
+    return __builtin_clz(hi);
+  else {
+    unsigned int low = n;
+    return (WORD_SIZE >> 1) + __builtin_clz(low);
+  }
+# else
+  return __builtin_clz(n);
+# endif
+#else
+  int c = 0, d = (WORD_SIZE >> 1);
+  while (d) {
+    if (n >> (c + d))
+      c += d;
+    d = d >> 1;
+  }
+  return WORD_SIZE - 1 - c;
+#endif
+}
+
+XFORM_NONGCING static int any_nonzero_digits(bigdig *na, intptr_t nl, int delta)
+/* if `delta`, then check only after that many bits in the most-significant
+   digit */
+{
+  if (delta) {
+    if (na[nl-1] & (((bigdig)1 << (WORD_SIZE - delta)) - 1))
+      return 1;
+    nl--;
+  }
+  
+  while (nl--) {
+    if (na[nl])
+      return 1;
+  }
+  return 0;
+}
+
 #define USE_FLOAT_BITS 53
 #define FP_TYPE double
 
@@ -1431,6 +1473,7 @@ static void bignum_add1_inplace(Scheme_Object **_stk_o)
 #define FP_TYPE_MULT(x, y) ((x)*(y))
 #define FP_TYPE_PLUS(x, y) ((x)+(y))
 #define FP_TYPE_DIV(x, y) ((x)/(y))
+#define FP_TYPE_POW(x, y) pow(x, y)
 #define FP_TYPE_FROM_INT(x) ((FP_TYPE)(x))
 #define FP_TYPE_GREATER_OR_EQV(x, y) ((x)>=(y))
 #define FP_TYPE_MINUS(x, y) ((x)-(y))
@@ -1453,6 +1496,7 @@ static void bignum_add1_inplace(Scheme_Object **_stk_o)
 #define FP_TYPE_MULT(x, y) ((x)*(y))
 #define FP_TYPE_PLUS(x, y) ((x)+(y))
 #define FP_TYPE_DIV(x, y) ((x)/(y))
+#define FP_TYPE_POW(x, y) pow(x, y)
 #define FP_TYPE_FROM_INT(x) ((FP_TYPE)(x))
 #define FP_TYPE_GREATER_OR_EQV(x, y) ((x)>=(y))
 #define FP_TYPE_MINUS(x, y) ((x)-(y))
@@ -1475,6 +1519,7 @@ static void bignum_add1_inplace(Scheme_Object **_stk_o)
 # define FP_TYPE_MULT(x, y) long_double_mult(x, y)
 # define FP_TYPE_DIV(x, y) long_double_div(x, y)
 # define FP_TYPE_PLUS(x, y) long_double_plus(x, y)
+# define FP_TYPE_POW(x, y) long_double_pow(x, y)
 # define FP_TYPE_FROM_INT(x) long_double_from_int(x)
 # define FP_TYPE_GREATER_OR_EQV(x, y) long_double_greater_or_eqv(x, y)
 # define FP_TYPE_MINUS(x, y) long_double_minus(x, y)