From 54f194f285114ad046fb67fd3027053f9544080f Mon Sep 17 00:00:00 2001
From: Mark H Weaver <mhw@netris.org>
Date: Fri, 4 Mar 2011 15:28:51 -0500
Subject: [PATCH 5/7] Simplify and improve scm_i_big2dbl, and add scm_i_big2dbl_2exp

* libguile/numbers.c (scm_i_big2dbl_2exp): New internal function which
  converts a bignum into a normalized significand and an exponent,
  similar to frexp, or mpz_get_d_2exp but with correct rounding.  This
  can be used on bignums that are too large or small to be represented
  by a double.

  (scm_i_big2dbl): Simplify and improve scm_i_big2dbl, which is now
  implemented in terms of scm_i_big2dbl_2exp.  Ties now go to even.
  Previously they were rounded away from zero.
---
 libguile/numbers.c |   94 ++++++++++++++++-----------------------------------
 1 files changed, 30 insertions(+), 64 deletions(-)

diff --git a/libguile/numbers.c b/libguile/numbers.c
index 1274306..1773715 100644
--- a/libguile/numbers.c
+++ b/libguile/numbers.c
@@ -331,76 +331,42 @@ 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.
-
-   The approach is to use mpz_get_d to pick out the high
-   double_precision 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.
-
-   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 double_precision 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)
-{
-  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 */
-    mpz_t  tmp;
-    if (bits > double_precision)
-      {
-        size_t  shift = bits - double_precision;
-        mpz_init2 (tmp, double_precision);
-        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
+static SCM round_right_shift_exact_integer (SCM n, long count);
+
+/* 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. */
+static double
+scm_i_big2dbl_2exp (SCM b, long *expon_p)
+{
+  size_t bits = mpz_sizeinbase (SCM_I_BIG_MPZ (b), 2);
+  size_t shift = 0;
+  double signif;
+  long expon;
 
   if (bits > double_precision)
     {
-      unsigned long  pos = bits - double_precision - 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)))
-        {
-          result += ldexp ((double) mpz_sgn (SCM_I_BIG_MPZ (b)), pos + 1);
-        }
+      shift = bits - double_precision;
+      b = round_right_shift_exact_integer (b, shift);
     }
-
+  signif = mpz_get_d_2exp (&expon, SCM_I_BIG_MPZ (b));
   scm_remember_upto_here_1 (b);
-  return result;
+  *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
-- 
1.7.2.5