From 04c4538f44a9222fa6bf1b7cdfa366982419a7d3 Mon Sep 17 00:00:00 2001 From: Matthew Flatt Date: Sun, 27 Apr 2014 15:17:57 -0600 Subject: [PATCH] speed up `inexact->exact` Parse IEEE floating-point numbers directly. Thanks to Neil T.! --- .../racket-test/tests/racket/number.rktl | 33 ++- racket/src/racket/sconfig.h | 3 + racket/src/racket/src/numstr.c | 205 +++++++++++------- racket/src/racket/src/ratfloat.inc | 73 ++++++- racket/src/racket/src/rational.c | 18 +- racket/src/racket/src/schpriv.h | 2 + 6 files changed, 236 insertions(+), 98 deletions(-) diff --git a/pkgs/racket-pkgs/racket-test/tests/racket/number.rktl b/pkgs/racket-pkgs/racket-test/tests/racket/number.rktl index d295a3f453..601d781166 100644 --- a/pkgs/racket-pkgs/racket-test/tests/racket/number.rktl +++ b/pkgs/racket-pkgs/racket-test/tests/racket/number.rktl @@ -3082,23 +3082,25 @@ (require racket/extflonum) -(define (check start end exact-> >=?) +(define (check start end exact-> ->exact >=?) (define delta (/ (- end start) 300)) (for/fold ([prev (exact-> start)]) ([i (in-range start (+ end delta) delta)]) (define next (exact-> i)) (test #t >=? next prev) + (test next exact-> (->exact next)) next) (for/fold ([prev (exact-> start)]) ([i (in-range start (+ end delta) delta)]) (define next (exact-> (- i))) (test #t >=? prev next) + (test next exact-> (->exact next)) next) (void)) -(check #e100000000000000.0 #e100000000000000.1 exact->inexact >=) -(check #e100000000000000.0 #e100000000000000.1 real->double-flonum >=) -(check #e1000000.0 #e1000000.1 real->single-flonum >=) -(when extflonum-available? - (check #e1000000000000000000.0 #e1000000000000000000.1 real->extfl extfl>=)) +(check #e100000000000000.0 #e100000000000000.1 exact->inexact inexact->exact >=) +(check #e100000000000000.0 #e100000000000000.1 real->double-flonum inexact->exact >=) +(check #e1000000.0 #e1000000.1 real->single-flonum inexact->exact >=) +(when (extflonum-available?) + (check #e1000000000000000000.0 #e1000000000000000000.1 real->extfl extfl->exact extfl>=)) ;; Sanity check (test 0.14285714285714285 real->double-flonum 1/7) @@ -3111,8 +3113,6 @@ (when (collection-file-path "base.rkt" "math" #:fail (lambda (x) #f)) (eval '(begin - (test #t string? "Randomized testing of rational->flonum") - (require math/base math/flonum) @@ -3123,12 +3123,27 @@ (* (if (< (random) 0.5) -1 1) (/ (random-bits (+ 1 (random 8192))) d))])) + (test #t string? "Randomized testing of rational->flonum") (for ([_ (in-range 10000)]) (define ry (random-rational)) (define y (real->double-flonum ry)) ; this generates rounding errors (define e (flulp-error y ry)) (unless (<= e 0.5) - (test #t (lambda (e y ry) (<= e 0.5)) e y ry)))))) + (test #t (lambda (e y ry) (<= e 0.5)) e y ry)) + (unless (= e (exact->inexact (inexact->exact e))) + (test e exact->inexact (inexact->exact e)))) + + (define -max (flonum->ordinal -max.0)) + (define +inf (flonum->ordinal +inf.0)) + (define (random-flonum) + (ordinal->flonum (random-integer -max +inf))) + + (test #t string? "Randomized testing of inexact->exact") + (define xs (list* -max.0 -1.0 -min.0 -0.0 0.0 +min.0 1.0 +max.0 + (build-list 100000 (λ (_) (random-flonum))))) + (for ([x (in-list xs)]) + (unless (= x (exact->inexact (inexact->exact x))) + (test x exact->inexact (inexact->exact x))))))) ;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; diff --git a/racket/src/racket/sconfig.h b/racket/src/racket/sconfig.h index 8ef4db54c1..8ab7e7a42c 100644 --- a/racket/src/racket/sconfig.h +++ b/racket/src/racket/sconfig.h @@ -1354,6 +1354,9 @@ /* Inexact Arithmetic */ /**********************/ + /* FLOATING_POINT_IS_NOT_IEEE disables inexact->exact conversion via + parsing of IEEE-format bits. */ + /* USE_SINGLE_FLOATS turns on support for single-precision floating point numbers. Otherwise, floating point numbers are always represented in double-precision. */ diff --git a/racket/src/racket/src/numstr.c b/racket/src/racket/src/numstr.c index 6ec8ea4ff9..b733697358 100644 --- a/racket/src/racket/src/numstr.c +++ b/racket/src/racket/src/numstr.c @@ -2036,6 +2036,130 @@ int scheme_check_long_double(const char *where, long_double d, const char *dest) /* native representations */ /*========================================================================*/ +Scheme_Object *scheme_bytes_to_integer(char *str, int slen, int sgned, int rshft, int mask) +{ + switch(slen) { + case 2: + if (sgned) { + short val; + memcpy(&val, str, sizeof(short)); + return scheme_make_integer(val); + } else { + unsigned short val; + memcpy(&val, str, sizeof(unsigned short)); + val >>= rshft; + if (mask < 16) { val &= (((unsigned short)1 << mask) - 1); } + return scheme_make_integer(val); + } + break; + case 4: + if (sgned) { + int val; + memcpy(&val, str, sizeof(int)); + return scheme_make_integer_value(val); + } else { + unsigned int val; + memcpy(&val, str, sizeof(unsigned int)); + val >>= rshft; + if (mask < 32) { val &= (((unsigned int)1 << mask) - 1); } + return scheme_make_integer_value_from_unsigned(val); + } + break; + default: +#ifdef SIXTY_FOUR_BIT_INTEGERS + if (sgned) { + intptr_t val; + memcpy(&val, str, sizeof(intptr_t)); + return scheme_make_integer_value(val); + } + else { + uintptr_t val; + memcpy(&val, str, sizeof(uintptr_t)); + val >>= rshft; + if (mask < 64) { val &= (((uintptr_t)1 << mask) - 1); } + return scheme_make_integer_value_from_unsigned(val); + } + break; +#else +# ifndef NO_LONG_LONG_TYPE + { + if (sgned) { + mzlonglong lv; + memcpy(&lv, str, sizeof(mzlonglong)); + return scheme_make_integer_value_from_long_long(lv); + } else { + umzlonglong lv; + memcpy(&lv, str, sizeof(umzlonglong)); + lv >>= rshft; + if (mask < 64) { lv &= (((umzlonglong)1 << mask) - 1); } + return scheme_make_integer_value_from_unsigned_long_long(lv); + } + break; + } +# else + { + Scheme_Object *h, *l, *a[2]; + unsigned int val; + +# if MZ_IS_BIG_ENDIAN + /* make little-endian at int level: */ + { + int v; + v = ((int *)str)[0]; + buf[0] = ((int *)str)[1]; + buf[1] = v; + str = (char *)buf; + } +# endif + + if (rshft >= 32) { + + } + + if (sgned) + h = scheme_make_integer_value(((int *)str)[1]); + else { + memcpy(&val, str + sizeof(unsigned int), sizeof(unsigned int)); + if (rshft >= 32) { + rshft -= 32; + val >>= rshft; + if (mask < 64) { val &= (((umzlonglong)1 << mask) - 1); } + return scheme_make_integer_value_from_unsigned(val); + } else { + h = scheme_make_integer_value_from_unsigned(val); + } + } + + memcpy(&val, str, sizeof(unsigned int)); + val >>= rshft; + l = scheme_make_integer_value_from_unsigned(val); + + a[0] = h; + a[1] = scheme_make_integer(32-rshft); + h = scheme_bitwise_shift(2, a); + + l = scheme_bin_plus(h, l); + + if (mask < 64) { + a[0] = scheme_make_integer(1); + a[1] = scheme_make_integer(mask); + h = scheme_bitwise_shift(2, a); + h = scheme_bin_minus(h, scheme_make_integer(1)); + a[0] = h; + a[1] = l; + l = scheme_bitwise_and(2, a); + } + + return l; + } +# endif +#endif + break; + } + + ESCAPED_BEFORE_HERE; +} + static Scheme_Object *bytes_to_integer (int argc, Scheme_Object *argv[]) { intptr_t strlen, slen; @@ -2086,86 +2210,7 @@ static Scheme_Object *bytes_to_integer (int argc, Scheme_Object *argv[]) str = (char *)buf; } - switch(slen) { - case 2: - if (sgned) { - short val; - memcpy(&val, str, sizeof(short)); - return scheme_make_integer(val); - } - else { - unsigned short val; - memcpy(&val, str, sizeof(unsigned short)); - return scheme_make_integer(val); - } - break; - case 4: - if (sgned) { - int val; - memcpy(&val, str, sizeof(int)); - return scheme_make_integer_value(val); - } - else { - unsigned int val; - memcpy(&val, str, sizeof(unsigned int)); - return scheme_make_integer_value_from_unsigned(val); - } - break; - default: -#ifdef SIXTY_FOUR_BIT_INTEGERS - if (sgned) { - intptr_t val; - memcpy(&val, str, sizeof(intptr_t)); - return scheme_make_integer_value(val); - } - else { - uintptr_t val; - memcpy(&val, str, sizeof(uintptr_t)); - return scheme_make_integer_value_from_unsigned(val); - } - break; -#else -# ifndef NO_LONG_LONG_TYPE - { - mzlonglong lv; - memcpy(&lv, str, sizeof(mzlonglong)); - if (sgned) - return scheme_make_integer_value_from_long_long(lv); - else - return scheme_make_integer_value_from_unsigned_long_long((umzlonglong)lv); - break; - } -# else - { - Scheme_Object *h, *l, *a[2]; - -# if MZ_IS_BIG_ENDIAN - /* make little-endian at int level: */ - { - int v; - v = ((int *)str)[0]; - buf[0] = ((int *)str)[1]; - buf[1] = v; - str = (char *)buf; - } -# endif - - if (sgned) - h = scheme_make_integer_value(((int *)str)[1]); - else - h = scheme_make_integer_value_from_unsigned(((unsigned int *)str)[1]); - l = scheme_make_integer_value_from_unsigned(((unsigned int *)str)[0]); - a[0] = h; - a[1] = scheme_make_integer(32); - h = scheme_bitwise_shift(2, a); - return scheme_bin_plus(h, l); - } -# endif -#endif - break; - } - - /* throw an error here */ + return scheme_bytes_to_integer(str, slen, sgned, 0, slen*8); } #define MZ_U8HI 0 diff --git a/racket/src/racket/src/ratfloat.inc b/racket/src/racket/src/ratfloat.inc index 611848c211..cf486bea7c 100644 --- a/racket/src/racket/src/ratfloat.inc +++ b/racket/src/racket/src/ratfloat.inc @@ -151,17 +151,75 @@ FP_TYPE SCHEME_RATIONAL_TO_FLOAT(const Scheme_Object *o) Scheme_Object *SCHEME_RATIONAL_FROM_FLOAT(FP_TYPE d) { +#ifdef DECODE_IEEE_FLOATING_POINT + Scheme_Object *a[2], *r, *m; + intptr_t s, e; + + SCHEME_CHECK_FLOAT("inexact->exact", d, "exact"); + +# ifdef FP_FITS_INT_TYPE + { + FP_FITS_INT_TYPE k; + memcpy(&k, &d, sizeof(FP_TYPE)); + s = (k >> (FLOAT_M_BITS + FLOAT_E_BITS)) & 0x1; + e = (k >> FLOAT_M_BITS) & (((FP_FITS_INT_TYPE)1 << FLOAT_E_BITS) - 1); + m = scheme_make_integer(k & (((FP_FITS_INT_TYPE)1 << FLOAT_M_BITS) - 1)); + } +# else + { + char b[sizeof(FP_TYPE)]; + memcpy(b, &d, sizeof(FP_TYPE)); + m = scheme_bytes_to_integer(b, sizeof(FP_TYPE), 0, (FLOAT_M_BITS + FLOAT_E_BITS), 1); + s = SCHEME_INT_VAL(m); + m = scheme_bytes_to_integer(b, sizeof(FP_TYPE), 0, FLOAT_M_BITS, FLOAT_E_BITS); + e = SCHEME_INT_VAL(m); + m = scheme_bytes_to_integer(b, sizeof(FP_TYPE), 0, 0, FLOAT_M_BITS); + } +# endif + + if (e == 0) { + /* Subnormal numbers all have the same exponent, 2^FLOAT_E_MIN */ + a[0] = scheme_make_integer(1); + a[1] = scheme_make_integer(-FLOAT_E_MIN); + r = scheme_bin_div(m, scheme_bitwise_shift(2, a)); + } else { + /* Adjust exponent for bias and mantissa size; add implicit one bit to mantissa */ + e -= (((1 << (FLOAT_E_BITS-1))-1) + FLOAT_M_BITS); +# ifdef FP_FITS_INT_TYPE + m = scheme_make_integer(SCHEME_INT_VAL(m) | ((FP_FITS_INT_TYPE)1 << FLOAT_M_BITS)); +# else + a[0] = scheme_make_integer(1); + a[1] = scheme_make_integer(FLOAT_M_BITS); + m = scheme_bin_plus(m, scheme_bitwise_shift(2, a)); +#endif + /* Compute m * 2^e */ + if (e < 0) { + a[0] = scheme_make_integer(1); + a[1] = scheme_make_integer(-e); + r = scheme_bin_div(m, scheme_bitwise_shift(2, a)); + } else { + a[0] = m; + a[1] = scheme_make_integer(e); + r = scheme_bitwise_shift(2, a); + } + } + + if (s) + return scheme_bin_minus(scheme_make_integer(0), r); + else + return r; +#else FP_DOUBLE_TYPE frac, i; int count, exponent, is_neg; Scheme_Object *int_part, *frac_part, *frac_num, *frac_denom, *two, *result; -#ifdef COMPUTE_NEG_INEXACT_TO_EXACT_AS_POS +# ifdef COMPUTE_NEG_INEXACT_TO_EXACT_AS_POS int negate; if (d < FP_ZEROx) { d = -d; negate = 1; } else negate = 0; -#endif +# endif SCHEME_CHECK_FLOAT("inexact->exact", d, "exact"); @@ -173,10 +231,10 @@ Scheme_Object *SCHEME_RATIONAL_FROM_FLOAT(FP_TYPE d) int_part = SCHEME_BIGNUM_FROM_FLOAT(i); if (FP_EQV(frac, FP_ZEROx)) { -#ifdef COMPUTE_NEG_INEXACT_TO_EXACT_AS_POS +# ifdef COMPUTE_NEG_INEXACT_TO_EXACT_AS_POS if (negate) return scheme_bin_minus(scheme_make_integer(0), int_part); -#endif +# endif return int_part; } @@ -202,12 +260,13 @@ Scheme_Object *SCHEME_RATIONAL_FROM_FLOAT(FP_TYPE d) result = scheme_bin_plus(int_part, frac_part); -#ifdef COMPUTE_NEG_INEXACT_TO_EXACT_AS_POS +# ifdef COMPUTE_NEG_INEXACT_TO_EXACT_AS_POS if (negate) return scheme_bin_minus(scheme_make_integer(0), result); -#endif +# endif return result; +#endif } #undef FP_TYPE @@ -221,6 +280,7 @@ Scheme_Object *SCHEME_RATIONAL_FROM_FLOAT(FP_TYPE d) #undef DO_FLOAT_DIV #undef FLOAT_E_MIN #undef FLOAT_M_BITS +#undef FLOAT_E_BITS #undef FP_ZEROx #undef FP_POWx #undef FP_MODFx @@ -235,3 +295,4 @@ Scheme_Object *SCHEME_RATIONAL_FROM_FLOAT(FP_TYPE d) #undef FP_LDEXP #undef FP_EQV #undef FP_IS_ZERO +#undef FP_FITS_INT_TYPE diff --git a/racket/src/racket/src/rational.c b/racket/src/racket/src/rational.c index 5002fff1fa..c6674f4032 100644 --- a/racket/src/racket/src/rational.c +++ b/racket/src/racket/src/rational.c @@ -515,6 +515,10 @@ Scheme_Object *scheme_rational_sqrt(const Scheme_Object *o) #endif } +#ifndef FLOATING_POINT_IS_NOT_IEEE +# define DECODE_IEEE_FLOATING_POINT +#endif + #define FP_TYPE double #define FP_MULT(x, y) x*y #define FP_DIV(x, y) x/y @@ -526,6 +530,7 @@ Scheme_Object *scheme_rational_sqrt(const Scheme_Object *o) #ifdef SIXTY_FOUR_BIT_INTEGERS # define FIXNUM_FITS_FP(x) (!(SCHEME_INT_VAL(x) & ~(((intptr_t)1 << (FLOAT_M_BITS-1)) - 1))) # define BIGNUM_FITS_FP(x) 0 +# define FP_FITS_INT_TYPE uintptr_t #else # define FIXNUM_FITS_FP(x) 1 # define BIGNUM_FITS_FP(x) (scheme_integer_length(x) <= (FLOAT_M_BITS-1)) @@ -536,8 +541,9 @@ Scheme_Object *scheme_rational_sqrt(const Scheme_Object *o) #define SCHEME_CHECK_FLOAT scheme_check_double #define SCHEME_BIGNUM_FROM_FLOAT scheme_bignum_from_double #define DO_FLOAT_DIV scheme__do_double_div -#define FLOAT_E_MIN -1074 +#define FLOAT_E_MIN (-1074) #define FLOAT_M_BITS 52 +#define FLOAT_E_BITS 11 #include "ratfloat.inc" #ifdef MZ_USE_SINGLE_FLOATS @@ -551,17 +557,22 @@ Scheme_Object *scheme_rational_sqrt(const Scheme_Object *o) #define FIXNUM_FITS_FP(x) (!(SCHEME_INT_VAL(x) & ~(((intptr_t)1 << (FLOAT_M_BITS-1)) - 1))) #define FP_IS_ZERO(x) x==0.0 #define BIGNUM_FITS_FP(x) 0 +#define FP_FITS_INT_TYPE unsigned int #define SCHEME_RATIONAL_TO_FLOAT scheme_rational_to_float #define SCHEME_RATIONAL_FROM_FLOAT scheme_rational_from_float #define SCHEME_BIGNUM_TO_FLOAT_INF_INFO scheme_bignum_to_float_inf_info #define SCHEME_CHECK_FLOAT scheme_check_float #define SCHEME_BIGNUM_FROM_FLOAT scheme_bignum_from_float #define DO_FLOAT_DIV scheme__do_float_div -#define FLOAT_E_MIN -127 +#define FLOAT_E_MIN (-127) #define FLOAT_M_BITS 23 +#define FLOAT_E_BITS 8 #include "ratfloat.inc" #endif +/* scheme_bytes_to_integer() can't real with sizeof(long_double): */ +#undef DECODE_IEEE_FLOATING_POINT + #ifdef MZ_LONG_DOUBLE # define FP_TYPE long_double # define FP_MULT(x, y) long_double_mult(x,y) @@ -578,8 +589,9 @@ Scheme_Object *scheme_rational_sqrt(const Scheme_Object *o) # define SCHEME_BIGNUM_TO_FLOAT_INF_INFO scheme_bignum_to_long_double_inf_info # define SCHEME_CHECK_FLOAT scheme_check_long_double # define SCHEME_BIGNUM_FROM_FLOAT scheme_bignum_from_long_double -# define FLOAT_E_MIN -16383 +# define FLOAT_E_MIN (-16383) # define FLOAT_M_BITS 64 +# define FLOAT_E_BITS 15 # define FP_ZEROx get_long_double_zero() # define FP_POWx pow # define FP_MODFx long_double_modf diff --git a/racket/src/racket/src/schpriv.h b/racket/src/racket/src/schpriv.h index cef4914240..be48d6a773 100644 --- a/racket/src/racket/src/schpriv.h +++ b/racket/src/racket/src/schpriv.h @@ -2259,6 +2259,8 @@ int scheme_nonneg_exact_p(Scheme_Object *n); Scheme_Object *scheme_floor(int argc, Scheme_Object *argv[]); +Scheme_Object *scheme_bytes_to_integer(char *str, int slen, int sgned, int rshft, int mask); + #ifdef TIME_TYPE_IS_UNSIGNED # define scheme_make_integer_value_from_time(t) scheme_make_integer_value_from_unsigned((uintptr_t)t) # define scheme_get_time_val(o, v) scheme_get_unsigned_int_val(o, v)