suzanne.soy/racket
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

speed up inexact->exact

Browse Source 
Parse IEEE floating-point numbers directly.

Thanks to Neil T.!
This commit is contained in:
Matthew Flatt 2014-04-27 15:17:57 -06:00
parent d682c940bd
commit 04c4538f44
6 changed files with 236 additions and 98 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

33
pkgs/racket-pkgs/racket-test/tests/racket/number.rktl
View File

@ -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)))))))
;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;

3
racket/src/racket/sconfig.h
View File

@ -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. */

205
racket/src/racket/src/numstr.c
View File

@ -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

73
racket/src/racket/src/ratfloat.inc
View File

@ -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

18
racket/src/racket/src/rational.c
View File

@ -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

2
racket/src/racket/src/schpriv.h
View File

@ -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)