Reimplement 'inexact->exact' to avoid mpq functions.

* libguile/numbers.c (scm_inexact_to_exact): Implement conversion of a
  double to an exact rational without using the mpq functions.

* test-suite/tests/numbers.test (dbl-mant-dig): Simplify initializer.
  (dbl-epsilon, dbl-min-exp): New variables.
  ("inexact->exact"): Add tests.  Fix broken "2.0**i to exact and back"
  test, and change it to "2.0**i to exact", to avoid use of
  'exact->inexact'.

diff --git a/libguile/numbers.c b/libguile/numbers.c
index fa55b4f..f0f7236 100644 (file)
--- a/libguile/numbers.c
+++ b/libguile/numbers.c
@@ -9109,22 +9109,35 @@ SCM_PRIMITIVE_GENERIC (scm_inexact_to_exact, "inexact->exact", 1, 0, 0,
 
       if (!SCM_LIKELY (DOUBLE_IS_FINITE (val)))
        SCM_OUT_OF_RANGE (1, z);
+      else if (val == 0.0)
+        return SCM_INUM0;
       else
        {
-         mpq_t frac;
-         SCM q;
-         
-         mpq_init (frac);
-         mpq_set_d (frac, val);
-         q = scm_i_make_ratio_already_reduced
-           (scm_i_mpz2num (mpq_numref (frac)),
-            scm_i_mpz2num (mpq_denref (frac)));
-
-         /* When scm_i_make_ratio throws, we leak the memory allocated
-            for frac...
-          */
-         mpq_clear (frac);
-         return q;
+          int expon;
+          SCM numerator;
+
+          numerator = scm_i_dbl2big (ldexp (frexp (val, &expon),
+                                            DBL_MANT_DIG));
+          expon -= DBL_MANT_DIG;
+          if (expon < 0)
+            {
+              int shift = mpz_scan1 (SCM_I_BIG_MPZ (numerator), 0);
+
+              if (shift > -expon)
+                shift = -expon;
+              mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (numerator),
+                               SCM_I_BIG_MPZ (numerator),
+                               shift);
+              expon += shift;
+            }
+          numerator = scm_i_normbig (numerator);
+          if (expon < 0)
+            return scm_i_make_ratio_already_reduced
+              (numerator, left_shift_exact_integer (SCM_INUM1, -expon));
+          else if (expon > 0)
+            return left_shift_exact_integer (numerator, expon);
+          else
+            return numerator;
        }
     }
 }

diff --git a/test-suite/tests/numbers.test b/test-suite/tests/numbers.test
index 6b4e08c..550dc50 100644 (file)
--- a/test-suite/tests/numbers.test
+++ b/test-suite/tests/numbers.test
@@ -46,15 +46,24 @@
 ;; the usual 53.
 ;;
 (define dbl-mant-dig
-  (let more ((i 1)
-            (d 2.0))
-    (if (> i 1024)
-       (error "Oops, cannot determine number of bits in mantissa of inexact"))
-    (let* ((sum  (+ 1.0 d))
-          (diff (- sum d)))
-      (if (= diff 1.0)
-         (more (1+ i) (* 2.0 d))
-         i))))
+  (do ((prec 0 (+ prec 1))
+       (eps 1.0 (/ eps 2.0)))
+      ((begin (when (> prec 1000000)
+                (error "Unable to determine dbl-mant-dig"))
+              (= 1.0 (+ 1.0 eps)))
+       prec)))
+
+(define dbl-epsilon
+  (expt 0.5 (- dbl-mant-dig 1)))
+
+(define dbl-min-exp
+  (do ((x 1.0 (/ x 2.0))
+       (y (+ 1.0 dbl-epsilon) (/ y 2.0))
+       (e 2 (- e 1)))
+      ((begin (when (< e -100000000)
+                (error "Unable to determine dbl-min-exp"))
+              (= x y))
+       e)))
 
 ;; like ash, but working on a flonum
 (define (ash-flo x n)
@@ -4251,6 +4260,13 @@
 ;;;
 
 (with-test-prefix "inexact->exact"
+
+  ;; Test "(inexact->exact f)", expect "want".
+  (define (test name f want)
+    (with-test-prefix name
+      (pass-if-equal "pos" want (inexact->exact f))
+      (pass-if-equal "neg" (- want) (inexact->exact (- f)))))
+
   (pass-if (documented? inexact->exact))
 
   (pass-if-exception "+inf" exception:out-of-range
@@ -4261,13 +4277,49 @@
   
   (pass-if-exception "nan" exception:out-of-range
     (inexact->exact +nan.0))
-  
-  (with-test-prefix "2.0**i to exact and back"
+
+  (test "0.0" 0.0 0)
+  (test "small even integer" 72.0 72)
+  (test "small odd integer"  73.0 73)
+
+  (test "largest inexact odd integer"
+        (- (expt 2.0 dbl-mant-dig) 1)
+        (- (expt 2   dbl-mant-dig) 1))
+
+  (test "largest inexact odd integer - 1"
+        (- (expt 2.0 dbl-mant-dig) 2)
+        (- (expt 2   dbl-mant-dig) 2))
+
+  (test "largest inexact odd integer + 3"
+        (+ (expt 2.0 dbl-mant-dig) 2)
+        (+ (expt 2   dbl-mant-dig) 2))
+
+  (test "largest inexact odd integer * 2^48"
+        (* (expt 2.0 48) (- (expt 2.0 dbl-mant-dig) 1))
+        (* (expt 2   48) (- (expt 2   dbl-mant-dig) 1)))
+
+  (test "largest inexact odd integer / 2^48"
+        (* (expt 0.5 48) (- (expt 2.0 dbl-mant-dig) 1))
+        (* (expt 1/2 48) (- (expt 2   dbl-mant-dig) 1)))
+
+  (test "smallest inexact"
+        (expt 2.0 (- dbl-min-exp dbl-mant-dig))
+        (expt 2   (- dbl-min-exp dbl-mant-dig)))
+
+  (test "smallest inexact * 2"
+        (* 2.0 (expt 2.0 (- dbl-min-exp dbl-mant-dig)))
+        (* 2   (expt 2   (- dbl-min-exp dbl-mant-dig))))
+
+  (test "smallest inexact * 3"
+        (* 3.0 (expt 2.0 (- dbl-min-exp dbl-mant-dig)))
+        (* 3   (expt 2   (- dbl-min-exp dbl-mant-dig))))
+
+  (with-test-prefix "2.0**i to exact"
     (do ((i 0   (1+ i))
-        (n 1.0 (* 2.0 n)))
+         (n 1   (* 2 n))
+        (f 1.0 (* 2.0 f)))
        ((> i 100))
-      (pass-if (list i n)
-       (= n (inexact->exact (exact->inexact n)))))))
+      (test (list i n) f n))))
 
 ;;;
 ;;; integer-expt