From 6e6991d9f7043aa70dc53c12afb97eeaf3cd63ba Mon Sep 17 00:00:00 2001
From: Mark H Weaver <mhw@netris.org>
Date: Mon, 4 Mar 2013 18:42:27 -0500
Subject: [PATCH 5/5] 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'.
---
 libguile/numbers.c            |   41 +++++++++++++--------
 test-suite/tests/numbers.test |   80 +++++++++++++++++++++++++++++++++--------
 2 files changed, 93 insertions(+), 28 deletions(-)

diff --git a/libguile/numbers.c b/libguile/numbers.c
index 0b13d69..845b37a 100644
--- a/libguile/numbers.c
+++ b/libguile/numbers.c
@@ -9085,22 +9085,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 58ff7b8..93caf04 100644
--- 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)
@@ -4241,6 +4250,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
@@ -4251,13 +4267,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
-- 
1.7.10.4