diff --git a/pkgs/racket-test-core/tests/racket/number.rktl b/pkgs/racket-test-core/tests/racket/number.rktl index 00b73fbbf6..a00c9724a6 100644 --- a/pkgs/racket-test-core/tests/racket/number.rktl +++ b/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) diff --git a/racket/collects/compiler/private/xform.rkt b/racket/collects/compiler/private/xform.rkt index 38052342fe..ce0b2d9559 100644 --- a/racket/collects/compiler/private/xform.rkt +++ b/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 diff --git a/racket/src/configure b/racket/src/configure index 21aecfa370..6e6360cb35 100755 --- a/racket/src/configure +++ b/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 diff --git a/racket/src/racket/configure.ac b/racket/src/racket/configure.ac index bff753b8ed..579b400c39 100644 --- a/racket/src/racket/configure.ac +++ b/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 diff --git a/racket/src/racket/mzconfig.h.in b/racket/src/racket/mzconfig.h.in index 8987178abf..c110bdd358 100644 --- a/racket/src/racket/mzconfig.h.in +++ b/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 diff --git a/racket/src/racket/src/bgnfloat.inc b/racket/src/racket/src/bgnfloat.inc index bd07924ae9..38bd0ec5d1 100644 --- a/racket/src/racket/src/bgnfloat.inc +++ b/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 diff --git a/racket/src/racket/src/bignum.c b/racket/src/racket/src/bignum.c index 51ec30b937..f660b2a62b 100644 --- a/racket/src/racket/src/bignum.c +++ b/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)