repairs to precision of exact->inexact et al.

Thanks to Neil T.!

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

@@ -3077,6 +3077,59 @@
   
   (test #t list? (filter n-digit-has-nth-root? (build-list 5000 (lambda (x) (+ x 1))))))
 
+;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
+;; exact->inexact precision (thanks to Neil Toronto)
+
+(require racket/extflonum)
+
+(define (check start end 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)
+    next)
+  (for/fold ([prev (exact-> start)]) ([i (in-range start (+ end delta) delta)])
+    (define next (exact-> (- i)))
+    (test #t >=? prev 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>=))
+
+;; Sanity check
+(test 0.14285714285714285 real->double-flonum 1/7)
+(test 1.2857142857142858 real->double-flonum 9/7)
+;; Cases that real->double-flonum used to get wrong
+(test -4882.526517254422 real->double-flonum -13737024017780747/2813507303900)
+(test -9.792844933246106e-14 real->double-flonum -1656/16910305547451097)
+
+;; Hack to use the "math" package when it's available:
+(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)
+
+      (define (random-rational)
+        (define d (random-bits (+ 1 (random 8192))))
+        (cond [(zero? d)  (random-rational)]
+              [else
+               (* (if (< (random) 0.5) -1 1)
+                  (/ (random-bits (+ 1 (random 8192))) d))]))
+
+      (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))))))
+
 ;; ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;
 
 (report-errs)

racket/src/racket/src/numarith.c

@@ -1103,21 +1103,39 @@ scheme_modulo(int argc, Scheme_Object *argv[])
   return rem_mod(argc, argv, "modulo", 0);
 }
 
+static Scheme_Object *
+do_quotient_remainder(const Scheme_Object *n1, const Scheme_Object *n2, Scheme_Object **_rem)
+{
+  Scheme_Object *rem = NULL, *quot, *a[2];
+
+  quot = do_bin_quotient("quotient/remainder", n1, n2, &rem);
+  if (!rem) {
+    a[0] = (Scheme_Object *)n1;
+    a[1] = (Scheme_Object *)n2;
+    rem = rem_mod(2, a, "remainder", 1);
+  }
+  *_rem = rem;
+
+  return quot;
+}
 
 Scheme_Object *
 quotient_remainder(int argc, Scheme_Object *argv[])
 {
-  Scheme_Object *rem = NULL, *quot, *a[2];
+  Scheme_Object *rem, *quot, *a[2];
 
-  quot = do_bin_quotient("quotient/remainder", argv[0], argv[1], &rem);
-  if (!rem) {
-    rem = rem_mod(argc, argv, "remainder", 1);
-  }
+  quot = do_quotient_remainder(argv[0], argv[1], &rem);
   a[0] = quot;
   a[1] = rem;
   return scheme_values(2, a);
 }
 
+Scheme_Object *scheme_bin_quotient_remainder(const Scheme_Object *n1, const Scheme_Object *n2,
+                                             Scheme_Object **_rem)
+{
+  return do_quotient_remainder(n1, n2, _rem);
+}
+
 /************************************************************************/
 /*                                Flfx                                  */
 /************************************************************************/

racket/src/racket/src/ratfloat.inc

@@ -27,26 +27,48 @@ FP_TYPE SCHEME_RATIONAL_TO_FLOAT(const Scheme_Object *o)
   intptr_t ns, ds;
 
   if (SCHEME_INTP(r->num)) {
-  #ifdef CONVERT_INT_TO_FLOAT
+    if (FIXNUM_FITS_FP(r->num)) {
+#ifdef CONVERT_INT_TO_FLOAT
       n = CONVERT_INT_TO_FLOAT(SCHEME_INT_VAL(r->num));
-  #else
+#else
       n = FP_TYPE_FROM_INT(SCHEME_INT_VAL(r->num));
-  #endif
+#endif
       ns = 0;
-  } else
+    } else {
+      n = FP_ZEROx;
+      ns = 1;
+    }
+  } else {
+    if (BIGNUM_FITS_FP(r->num)) {
       n = SCHEME_BIGNUM_TO_FLOAT_INF_INFO(r->num, 0, &ns);
+    } else {
+      n = FP_ZEROx;
+      ns = 1;
+    }
+  }
 
   if (SCHEME_INTP(r->denom)) {
+    if (FIXNUM_FITS_FP(r->denom)) {
       d = FP_TYPE_FROM_INT(SCHEME_INT_VAL(r->denom));
       ds = 0;
-  } else
+    } else {
+      d = FP_ZEROx;
+      ds = 1;
+    }
+  } else {
+    if (BIGNUM_FITS_FP(r->denom)) {
       d = SCHEME_BIGNUM_TO_FLOAT_INF_INFO(r->denom, 0, &ds);
+    } else {
+      d = FP_ZEROx;
+      ds = 1;
+    }
+  }
 
   if (ns || ds) {
     /* Quick path doesn't necessarily work. The more general
        way is adpated from Gambit-C 4.1. */
     intptr_t nl, dl, p, shift;
-    Scheme_Object *a[2], *n, *d;
+    Scheme_Object *a[2], *n, *d, *rem;
     FP_TYPE res;
 
     a[0] = r->num;
@@ -87,9 +109,25 @@ FP_TYPE SCHEME_RATIONAL_TO_FLOAT(const Scheme_Object *o)
     a[1] = scheme_make_integer(shift);
     n = scheme_bitwise_shift(2, a);
 
-    n = scheme_bin_div(n, d);
-    if (SCHEME_RATIONALP(n))
-      n = scheme_rational_round(n);
+    /* Rounded divide: */
+    n = scheme_bin_quotient_remainder(n, d, &rem);
+    a[0] = d;
+    a[1] = scheme_make_integer(-1);
+    d = scheme_bitwise_shift(2, a);
+    if (!scheme_bin_lt(rem, d)) {
+      if (scheme_bin_gt(rem, d)) {
+        n = scheme_bin_plus(n, scheme_make_integer(1));
+      } else {
+        /* Round to even: */
+        a[0] = d;
+        if (!scheme_odd_p(1, a)) {
+          a[0] = n;
+          if (!scheme_even_p(1, a)) {
+            n = scheme_bin_plus(n, scheme_make_integer(1));
+          }
+        }
+      }
+    }
 
     if (SCHEME_INTP(n))
       res = FP_TYPE_FROM_INT(SCHEME_INT_VAL(n));
@@ -178,6 +216,8 @@ Scheme_Object *SCHEME_RATIONAL_FROM_FLOAT(FP_TYPE d)
 #undef SCHEME_BIGNUM_TO_FLOAT_INF_INFO
 #undef SCHEME_BIGNUM_FROM_FLOAT
 #undef SCHEME_CHECK_FLOAT
+#undef FIXNUM_FITS_FP
+#undef BIGNUM_FITS_FP
 #undef DO_FLOAT_DIV 
 #undef FLOAT_E_MIN
 #undef FLOAT_M_BITS

racket/src/racket/src/rational.c

@@ -523,6 +523,13 @@ Scheme_Object *scheme_rational_sqrt(const Scheme_Object *o)
 #define FP_LESS(x, y) x<y
 #define FP_IS_ZERO(x) x==0.0
 #define FP_TYPE_FROM_INT(x) (FP_TYPE)x
+#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
+#else
+# define FIXNUM_FITS_FP(x) 1
+# define BIGNUM_FITS_FP(x) (scheme_integer_length(x) <= (FLOAT_M_BITS-1))
+#endif
 #define SCHEME_RATIONAL_TO_FLOAT scheme_rational_to_double
 #define SCHEME_RATIONAL_FROM_FLOAT scheme_rational_from_double
 #define SCHEME_BIGNUM_TO_FLOAT_INF_INFO scheme_bignum_to_double_inf_info
@@ -541,7 +548,9 @@ Scheme_Object *scheme_rational_sqrt(const Scheme_Object *o)
 #define FP_EQV(x,y) x==y
 #define FP_LESS(x, y) x<y
 #define FP_TYPE_FROM_INT(x) (FP_TYPE)x
+#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 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
@@ -561,6 +570,8 @@ Scheme_Object *scheme_rational_sqrt(const Scheme_Object *o)
 # define FP_EQV(x,y) long_double_eqv(x,y)
 # define FP_LESS(x, y) long_double_less(x,y)
 # define FP_TYPE_FROM_INT(x) long_double_from_int(x)
+# define FIXNUM_FITS_FP(x) 1
+# define BIGNUM_FITS_FP(x) (scheme_integer_length(x) <= (FLOAT_M_BITS-1))
 # define FP_IS_ZERO(x) long_double_is_zero(x)
 # define SCHEME_RATIONAL_TO_FLOAT scheme_rational_to_long_double
 # define SCHEME_RATIONAL_FROM_FLOAT scheme_rational_from_long_double

racket/src/racket/src/schpriv.h

@@ -2226,6 +2226,8 @@ int scheme_bin_gt(const Scheme_Object *n1, const Scheme_Object *n2);
 int scheme_bin_gt_eq(const Scheme_Object *n1, const Scheme_Object *n2);
 int scheme_bin_lt_eq(const Scheme_Object *n1, const Scheme_Object *n2);
 
+Scheme_Object *scheme_bin_quotient_remainder(const Scheme_Object *n1, const Scheme_Object *n2, Scheme_Object **_rem);
+
 Scheme_Object *scheme_sub1(int argc, Scheme_Object *argv[]);
 Scheme_Object *scheme_add1(int argc, Scheme_Object *argv[]);
 Scheme_Object *scheme_odd_p(int argc, Scheme_Object *argv[]);