Simplify and improve scm_i_big2dbl, and add scm_i_big2dbl_2exp

* libguile/numbers.c (scm_i_big2dbl_2exp): New static function.
  (scm_i_big2dbl): Reimplement in terms of 'scm_i_big2dbl_2exp',
  with proper rounding.

* test-suite/tests/numbers.test ("exact->inexact"): Add tests.

2 files changed, 80 insertions, 78 deletions

diff --git a/libguile/numbers.c b/libguile/numbers.c
index 3f2afdebb..81461792d 100644
--- a/libguile/numbers.c
+++ b/libguile/numbers.c
@@ -330,81 +330,52 @@ scm_i_dbl2num (double u)
     return scm_i_dbl2big (u);
 }
 
-/* scm_i_big2dbl() rounds to the closest representable double, in accordance
-   with R5RS exact->inexact.
+static SCM round_right_shift_exact_integer (SCM n, long count);
 
-   The approach is to use mpz_get_d to pick out the high DBL_MANT_DIG bits
-   (ie. truncate towards zero), then adjust to get the closest double by
-   examining the next lower bit and adding 1 (to the absolute value) if
-   necessary.
+/* scm_i_big2dbl_2exp() is like frexp for bignums: it converts the
+   bignum b into a normalized significand and exponent such that
+   b = significand * 2^exponent and 1/2 <= abs(significand) < 1.
+   The return value is the significand rounded to the closest
+   representable double, and the exponent is placed into *expon_p.
+   If b is zero, then the returned exponent and significand are both
+   zero. */
 
-   Bignums exactly half way between representable doubles are rounded to the
-   next higher absolute value (ie. away from zero).  This seems like an
-   adequate interpretation of R5RS "numerically closest", and it's easier
-   and faster than a full "nearest-even" style.
-
-   The bit test must be done on the absolute value of the mpz_t, which means
-   we need to use mpz_getlimbn.  mpz_tstbit is not right, it treats
-   negatives as twos complement.
-
-   In GMP before 4.2, mpz_get_d rounding was unspecified.  It ended up
-   following the hardware rounding mode, but applied to the absolute
-   value of the mpz_t operand.  This is not what we want so we put the
-   high DBL_MANT_DIG bits into a temporary.  Starting with GMP 4.2
-   (released in March 2006) mpz_get_d now always truncates towards zero.
-
-   ENHANCE-ME: The temporary init+clear to force the rounding in GMP
-   before 4.2 is a slowdown.  It'd be faster to pick out the relevant
-   high bits with mpz_getlimbn.  */
-
-double
-scm_i_big2dbl (SCM b)
+static double
+scm_i_big2dbl_2exp (SCM b, long *expon_p)
 {
-  double result;
-  size_t bits;
-
-  bits = mpz_sizeinbase (SCM_I_BIG_MPZ (b), 2);
-
-#if 1
-  {
-    /* For GMP earlier than 4.2, force truncation towards zero */
-
-    /* FIXME: DBL_MANT_DIG is the number of base-`FLT_RADIX' digits,
-       _not_ the number of bits, so this code will break badly on a
-       system with non-binary doubles.  */
-
-    mpz_t  tmp;
-    if (bits > DBL_MANT_DIG)
-      {
-        size_t  shift = bits - DBL_MANT_DIG;
-        mpz_init2 (tmp, DBL_MANT_DIG);
-        mpz_tdiv_q_2exp (tmp, SCM_I_BIG_MPZ (b), shift);
-        result = ldexp (mpz_get_d (tmp), shift);
-        mpz_clear (tmp);
-      }
-    else
-      {
-        result = mpz_get_d (SCM_I_BIG_MPZ (b));
-      }
-  }
-#else
-  /* GMP 4.2 or later */
-  result = mpz_get_d (SCM_I_BIG_MPZ (b));
-#endif
+  size_t bits = mpz_sizeinbase (SCM_I_BIG_MPZ (b), 2);
+  size_t shift = 0;
 
   if (bits > DBL_MANT_DIG)
     {
-      unsigned long  pos = bits - DBL_MANT_DIG - 1;
-      /* test bit number "pos" in absolute value */
-      if (mpz_getlimbn (SCM_I_BIG_MPZ (b), pos / GMP_NUMB_BITS)
-          & ((mp_limb_t) 1 << (pos % GMP_NUMB_BITS)))
+      shift = bits - DBL_MANT_DIG;
+      b = round_right_shift_exact_integer (b, shift);
+      if (SCM_I_INUMP (b))
         {
-          result += ldexp ((double) mpz_sgn (SCM_I_BIG_MPZ (b)), pos + 1);
+          int expon;
+          double signif = frexp (SCM_I_INUM (b), &expon);
+          *expon_p = expon + shift;
+          return signif;
         }
     }
 
-  scm_remember_upto_here_1 (b);
-  return result;
+  {
+    long expon;
+    double signif = mpz_get_d_2exp (&expon, SCM_I_BIG_MPZ (b));
+    scm_remember_upto_here_1 (b);
+    *expon_p = expon + shift;
+    return signif;
+  }
+}
+
+/* scm_i_big2dbl() rounds to the closest representable double,
+   in accordance with R5RS exact->inexact.  */
+double
+scm_i_big2dbl (SCM b)
+{
+  long expon;
+  double signif = scm_i_big2dbl_2exp (b, &expon);
+  return ldexp (signif, expon);
 }
 
 SCM
diff --git a/test-suite/tests/numbers.test b/test-suite/tests/numbers.test
index bb1424853..6b4e08c3a 100644
--- a/test-suite/tests/numbers.test
+++ b/test-suite/tests/numbers.test
@@ -3858,21 +3858,17 @@
 ;;;
 
 (with-test-prefix "exact->inexact"
-  
+
+  ;; Test "(exact->inexact n)", expect "want".
+  (define (test name n want)
+    (with-test-prefix name
+      (pass-if-equal "pos" want (exact->inexact n))
+      (pass-if-equal "neg" (- want) (exact->inexact (- n)))))
+
   ;; Test "(exact->inexact n)", expect "want".
   ;; "i" is a index, for diagnostic purposes.
   (define (try-i i n want)
-    (with-test-prefix (list i n want)
-      (with-test-prefix "pos"
-	(let ((got (exact->inexact n)))
-	  (pass-if "inexact?" (inexact? got))
-	  (pass-if (list "=" got) (= want got))))
-      (set! n    (- n))
-      (set! want (- want))
-      (with-test-prefix "neg"
-	(let ((got (exact->inexact n)))
-	  (pass-if "inexact?" (inexact? got))
-	  (pass-if (list "=" got) (= want got))))))
+    (test (list i n want) n want))
   
   (with-test-prefix "2^i, no round"
     (do ((i    0   (1+ i))
@@ -3945,7 +3941,42 @@
   ;; convert the num and den to doubles, resulting in infs.
   (pass-if "frac big/big, exceeding double"
     (let ((big (ash 1 4096)))
-      (= 1.0 (exact->inexact (/ (1+ big) big))))))
+      (= 1.0 (exact->inexact (/ (1+ big) big)))))
+
+  (test "round up to odd"
+        ;; =====================================================
+        ;; 11111111111111111111111111111111111111111111111111000101 ->
+        ;; 11111111111111111111111111111111111111111111111111001000
+        (+ (expt 2   (+ dbl-mant-dig 3)) -64 #b000101)
+        (+ (expt 2.0 (+ dbl-mant-dig 3)) -64 #b001000))
+
+  (test "round down to odd"
+        ;; =====================================================
+        ;; 11111111111111111111111111111111111111111111111111001011 ->
+        ;; 11111111111111111111111111111111111111111111111111001000
+        (+ (expt 2   (+ dbl-mant-dig 3)) -64 #b001011)
+        (+ (expt 2.0 (+ dbl-mant-dig 3)) -64 #b001000))
+
+  (test "round tie up to even"
+        ;; =====================================================
+        ;; 11111111111111111111111111111111111111111111111111011100 ->
+        ;; 11111111111111111111111111111111111111111111111111100000
+        (+ (expt 2   (+ dbl-mant-dig 3)) -64 #b011100)
+        (+ (expt 2.0 (+ dbl-mant-dig 3)) -64 #b100000))
+
+  (test "round tie down to even"
+        ;; =====================================================
+        ;; 11111111111111111111111111111111111111111111111111000100 ->
+        ;; 11111111111111111111111111111111111111111111111111000000
+        (+ (expt 2   (+ dbl-mant-dig 3)) -64 #b000100)
+        (+ (expt 2.0 (+ dbl-mant-dig 3)) -64 #b000000))
+
+  (test "round tie up to next power of two"
+        ;;  =====================================================
+        ;;  11111111111111111111111111111111111111111111111111111100 ->
+        ;; 100000000000000000000000000000000000000000000000000000000
+        (+ (expt 2 (+ dbl-mant-dig 3)) -64 #b111100)
+        (expt 2.0 (+ dbl-mant-dig 3))))
 
 ;;;
 ;;; expt