guile-devel
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[PATCHES] Numeric improvements


From: Mark H Weaver
Subject: [PATCHES] Numeric improvements
Date: Tue, 05 Mar 2013 06:04:11 -0500

Hello all,

Here are ten patches to improve numerics.  Among other things, this
eliminates the known obstacles to linking with mini-gmp
<http://bugs.gnu.org/10519>, fixes several problems with our number
printer <http://bugs.gnu.org/13757>, adds 'round-ash', and speeds up
handling of exact rationals.

Still to be done: add more tests to ensure full coverage of all code
paths.  It's possible that there are some bugs lurking.  However, I
wanted to post it now to allow review and possible work on mini-gmp
integration.

Comments and suggestions welcome.

     Mark


>From cd9784ed33d78e6647a752123bf7be91d65b5c96 Mon Sep 17 00:00:00 2001
From: Mark H Weaver <address@hidden>
Date: Sun, 3 Mar 2013 04:34:17 -0500
Subject: [PATCH 01/10] Improve code in scm_gcd for inum/inum case

* libguile/numbers.c (scm_gcd): Improve implementation of inum/inum case
  to be more clear and efficient.
---
 libguile/numbers.c |   54 +++++++++++++++++++++++++++++-----------------------
 1 file changed, 30 insertions(+), 24 deletions(-)

diff --git a/libguile/numbers.c b/libguile/numbers.c
index 66c95db..9c28a79 100644
--- a/libguile/numbers.c
+++ b/libguile/numbers.c
@@ -3889,52 +3889,58 @@ SCM_PRIMITIVE_GENERIC (scm_i_gcd, "gcd", 0, 2, 1,
 SCM
 scm_gcd (SCM x, SCM y)
 {
-  if (SCM_UNBNDP (y))
+  if (SCM_UNLIKELY (SCM_UNBNDP (y)))
     return SCM_UNBNDP (x) ? SCM_INUM0 : scm_abs (x);
   
-  if (SCM_I_INUMP (x))
+  if (SCM_LIKELY (SCM_I_INUMP (x)))
     {
-      if (SCM_I_INUMP (y))
+      if (SCM_LIKELY (SCM_I_INUMP (y)))
         {
           scm_t_inum xx = SCM_I_INUM (x);
           scm_t_inum yy = SCM_I_INUM (y);
           scm_t_inum u = xx < 0 ? -xx : xx;
           scm_t_inum v = yy < 0 ? -yy : yy;
           scm_t_inum result;
-          if (xx == 0)
+          if (SCM_UNLIKELY (xx == 0))
            result = v;
-         else if (yy == 0)
+         else if (SCM_UNLIKELY (yy == 0))
            result = u;
          else
            {
-             scm_t_inum k = 1;
-             scm_t_inum t;
+             int k = 0;
              /* Determine a common factor 2^k */
-             while (!(1 & (u | v)))
+             while (((u | v) & 1) == 0)
                {
-                 k <<= 1;
+                 k++;
                  u >>= 1;
                  v >>= 1;
                }
              /* Now, any factor 2^n can be eliminated */
-             if (u & 1)
-               t = -v;
+             if ((u & 1) == 0)
+               while ((u & 1) == 0)
+                 u >>= 1;
              else
+               while ((v & 1) == 0)
+                 v >>= 1;
+             /* Both u and v are now odd.  Subtract the smaller one
+                from the larger one to produce an even number, remove
+                more factors of two, and repeat. */
+             while (u != v)
                {
-                 t = u;
-               b3:
-                 t = SCM_SRS (t, 1);
+                 if (u > v)
+                   {
+                     u -= v;
+                     while ((u & 1) == 0)
+                       u >>= 1;
+                   }
+                 else
+                   {
+                     v -= u;
+                     while ((v & 1) == 0)
+                       v >>= 1;
+                   }
                }
-             if (!(1 & t))
-               goto b3;
-             if (t > 0)
-               u = t;
-             else
-               v = -t;
-             t = u - v;
-             if (t != 0)
-               goto b3;
-             result = u * k;
+             result = u << k;
            }
           return (SCM_POSFIXABLE (result)
                  ? SCM_I_MAKINUM (result)
-- 
1.7.10.4

>From f6201616f7304979a31ab41814e2b297f74a3484 Mon Sep 17 00:00:00 2001
From: Mark H Weaver <address@hidden>
Date: Sun, 3 Mar 2013 04:34:50 -0500
Subject: [PATCH 02/10] Optimize and simplify fractions code

* libguile/numbers.c (scm_exact_integer_quotient): New internal static
  function that computes the quotient of two exact integers when the
  remainder is known in advance to be zero.  For large integers this can
  be implemented more efficiently than when the remainder is unknown.

  (scm_i_make_ratio_already_reduced): New internal static function that
  creates a ratio when the numerator and denominator are already known
  to be reduced to lowest terms (i.e. when their gcd is 1).  This can be
  used in several places to avoid unnecessary gcd computations.

  (scm_i_make_ratio): Rewrite in terms of
  scm_i_make_ratio_already_reduced.  Don't bother checking to see if the
  denominator divides the numerator evenly.  This is wasted effort in
  the common case.  Instead, compute the gcd, reduce to lowest terms
  (using scm_exact_integer_quotient), and let
  scm_i_make_ratio_already_reduced do the integer check (by checking for
  a unit denominator).

  (scm_integer_expt): Optimize fraction case by (a/b)^n ==> (a^n)/(b^n),
  to avoid needless effort to reduce intermediate products to lowest
  terms.  If a and b have no common factors, then a^n and b^n have no
  common factors.  Use scm_i_make_ratio_already_reduced to construct the
  final result, so that no gcd computations are needed to exponentiate a
  fraction.

  (scm_abs, scm_magnitude, scm_difference): Use
  scm_i_make_ratio_already_reduced to avoid wasteful gcd computations
  when negating fractions.

  (do_divide): Use scm_i_make_ratio_already_reduced to avoid wasteful
  gcd computations when taking the reciprocal of a fraction or integer.

* test-suite/tests/numbers.test (expt, integer-expt): Add tests.
---
 libguile/numbers.c            |  242 +++++++++++++++++++++++++----------------
 test-suite/tests/numbers.test |    6 +
 2 files changed, 157 insertions(+), 91 deletions(-)

diff --git a/libguile/numbers.c b/libguile/numbers.c
index 9c28a79..b5ba1a0 100644
--- a/libguile/numbers.c
+++ b/libguile/numbers.c
@@ -439,96 +439,58 @@ scm_i_mpz2num (mpz_t b)
 /* this is needed when we want scm_divide to make a float, not a ratio, even 
if passed two ints */
 static SCM scm_divide2real (SCM x, SCM y);
 
+/* scm_i_make_ratio_already_reduced makes a ratio more efficiently,
+   when the following conditions are known in advance:
+
+    1. numerator and denominator are exact integers
+    2. numerator and denominator are reduced to lowest terms: gcd(n,d) == 1
+*/
 static SCM
-scm_i_make_ratio (SCM numerator, SCM denominator)
-#define FUNC_NAME "make-ratio"
+scm_i_make_ratio_already_reduced (SCM numerator, SCM denominator)
 {
-  /* First make sure the arguments are proper.
-   */
-  if (SCM_I_INUMP (denominator))
+  /* Flip signs so that the denominator is positive. */
+  if (scm_is_false (scm_positive_p (denominator)))
     {
-      if (scm_is_eq (denominator, SCM_INUM0))
+      if (SCM_UNLIKELY (scm_is_eq (denominator, SCM_INUM0)))
        scm_num_overflow ("make-ratio");
-      if (scm_is_eq (denominator, SCM_INUM1))
-       return numerator;
-    }
-  else 
-    {
-      if (!(SCM_BIGP(denominator)))
-       SCM_WRONG_TYPE_ARG (2, denominator);
-    }
-  if (!SCM_I_INUMP (numerator) && !SCM_BIGP (numerator))
-    SCM_WRONG_TYPE_ARG (1, numerator);
-
-  /* Then flip signs so that the denominator is positive.
-   */
-  if (scm_is_true (scm_negative_p (denominator)))
-    {
-      numerator = scm_difference (numerator, SCM_UNDEFINED);
-      denominator = scm_difference (denominator, SCM_UNDEFINED);
-    }
-
-  /* Now consider for each of the four fixnum/bignum combinations
-     whether the rational number is really an integer.
-  */
-  if (SCM_I_INUMP (numerator))
-    {
-      scm_t_inum x = SCM_I_INUM (numerator);
-      if (scm_is_eq (numerator, SCM_INUM0))
-       return SCM_INUM0;
-      if (SCM_I_INUMP (denominator))
+      else
        {
-         scm_t_inum y;
-         y = SCM_I_INUM (denominator);
-         if (x == y)
-           return SCM_INUM1;
-         if ((x % y) == 0)
-           return SCM_I_MAKINUM (x / y);
+         numerator = scm_difference (numerator, SCM_UNDEFINED);
+         denominator = scm_difference (denominator, SCM_UNDEFINED);
        }
-      else
-        {
-          /* When x == SCM_MOST_NEGATIVE_FIXNUM we could have the negative
-             of that value for the denominator, as a bignum.  Apart from
-             that case, abs(bignum) > abs(inum) so inum/bignum is not an
-             integer.  */
-          if (x == SCM_MOST_NEGATIVE_FIXNUM
-              && mpz_cmp_ui (SCM_I_BIG_MPZ (denominator),
-                             - SCM_MOST_NEGATIVE_FIXNUM) == 0)
-           return SCM_I_MAKINUM(-1);
-        }
     }
-  else if (SCM_BIGP (numerator))
+
+  /* Check for the integer case */
+  if (scm_is_eq (denominator, SCM_INUM1))
+    return numerator;
+
+  return scm_double_cell (scm_tc16_fraction,
+                         SCM_UNPACK (numerator),
+                         SCM_UNPACK (denominator), 0);
+}
+
+static SCM scm_exact_integer_quotient (SCM x, SCM y);
+
+static SCM
+scm_i_make_ratio (SCM numerator, SCM denominator)
+#define FUNC_NAME "make-ratio"
+{
+  /* Make sure the arguments are proper */
+  if (!SCM_LIKELY (SCM_I_INUMP (numerator) || SCM_BIGP (numerator)))
+    SCM_WRONG_TYPE_ARG (1, numerator);
+  else if (!SCM_LIKELY (SCM_I_INUMP (denominator) || SCM_BIGP (denominator)))
+    SCM_WRONG_TYPE_ARG (2, denominator);
+  else
     {
-      if (SCM_I_INUMP (denominator))
-       {
-         scm_t_inum yy = SCM_I_INUM (denominator);
-         if (mpz_divisible_ui_p (SCM_I_BIG_MPZ (numerator), yy))
-           return scm_divide (numerator, denominator);
-       }
-      else
+      SCM the_gcd = scm_gcd (numerator, denominator);
+      if (!(scm_is_eq (the_gcd, SCM_INUM1)))
        {
-         if (scm_is_eq (numerator, denominator))
-           return SCM_INUM1;
-         if (mpz_divisible_p (SCM_I_BIG_MPZ (numerator),
-                              SCM_I_BIG_MPZ (denominator)))
-           return scm_divide(numerator, denominator);
+         /* Reduce to lowest terms */
+         numerator = scm_exact_integer_quotient (numerator, the_gcd);
+         denominator = scm_exact_integer_quotient (denominator, the_gcd);
        }
+      return scm_i_make_ratio_already_reduced (numerator, denominator);
     }
-
-  /* No, it's a proper fraction.
-   */
-  {
-    SCM divisor = scm_gcd (numerator, denominator);
-    if (!(scm_is_eq (divisor, SCM_INUM1)))
-      {
-       numerator = scm_divide (numerator, divisor);
-       denominator = scm_divide (denominator, divisor);
-      }
-      
-    return scm_double_cell (scm_tc16_fraction,
-                           SCM_UNPACK (numerator),
-                           SCM_UNPACK (denominator), 0);
-  }
 }
 #undef FUNC_NAME
 
@@ -820,8 +782,9 @@ SCM_PRIMITIVE_GENERIC (scm_abs, "abs", 1, 0, 0,
     {
       if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (x))))
        return x;
-      return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x), 
SCM_UNDEFINED),
-                            SCM_FRACTION_DENOMINATOR (x));
+      return scm_i_make_ratio_already_reduced
+       (scm_difference (SCM_FRACTION_NUMERATOR (x), SCM_UNDEFINED),
+        SCM_FRACTION_DENOMINATOR (x));
     }
   else
     SCM_WTA_DISPATCH_1 (g_scm_abs, x, 1, s_scm_abs);
@@ -889,6 +852,83 @@ SCM_PRIMITIVE_GENERIC (scm_modulo, "modulo", 2, 0, 0,
 }
 #undef FUNC_NAME
 
+/* scm_exact_integer_quotient returns the exact integer q such that
+   n = q*d, for exact integers n and d, where d is known in advance to
+   divide n evenly (with zero remainder). */
+static SCM
+scm_exact_integer_quotient (SCM n, SCM d)
+#define FUNC_NAME "exact-integer-quotient"
+{
+  if (SCM_LIKELY (SCM_I_INUMP (n)))
+    {
+      scm_t_inum nn = SCM_I_INUM (n);
+      if (SCM_LIKELY (SCM_I_INUMP (d)))
+       {
+         scm_t_inum dd = SCM_I_INUM (d);
+         if (SCM_UNLIKELY (dd == 0))
+           scm_num_overflow ("exact-integer-quotient");
+         else
+           {
+             scm_t_inum qq = nn / dd;
+             if (SCM_LIKELY (SCM_FIXABLE (qq)))
+               return SCM_I_MAKINUM (qq);
+             else
+               return scm_i_inum2big (qq);
+           }
+       }
+      else if (SCM_LIKELY (SCM_BIGP (d)))
+       {
+         /* n is an inum and d is a bignum.  Given that d is known to
+            divide n evenly, there are only two possibilities: n is 0,
+            or else n is fixnum-min and d is abs(fixnum-min). */
+         if (nn == 0)
+           return SCM_INUM0;
+         else
+           return SCM_I_MAKINUM (-1);
+       }
+      else
+       SCM_WRONG_TYPE_ARG (2, d);
+    }
+  else if (SCM_LIKELY (SCM_BIGP (n)))
+    {
+      if (SCM_LIKELY (SCM_I_INUMP (d)))
+       {
+         scm_t_inum dd = SCM_I_INUM (d);
+         if (SCM_UNLIKELY (dd == 0))
+           scm_num_overflow ("exact-integer-quotient");
+         else if (SCM_UNLIKELY (dd == 1))
+           return n;
+         else
+           {
+             SCM q = scm_i_mkbig ();
+             if (dd > 0)
+               mpz_divexact_ui (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (n), dd);
+             else
+               {
+                 mpz_divexact_ui (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (n), -dd);
+                 mpz_neg (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (q));
+               }
+             scm_remember_upto_here_1 (n);
+             return scm_i_normbig (q);
+           }
+       }
+      else if (SCM_LIKELY (SCM_BIGP (d)))
+       {
+         SCM q = scm_i_mkbig ();
+         mpz_divexact (SCM_I_BIG_MPZ (q),
+                       SCM_I_BIG_MPZ (n),
+                       SCM_I_BIG_MPZ (d));
+         scm_remember_upto_here_2 (n, d);
+         return scm_i_normbig (q);
+       }
+      else
+       SCM_WRONG_TYPE_ARG (2, d);
+    }
+  else
+    SCM_WRONG_TYPE_ARG (1, n);
+}
+#undef FUNC_NAME
+
 /* two_valued_wta_dispatch_2 is a version of SCM_WTA_DISPATCH_2 for
    two-valued functions.  It is called from primitive generics that take
    two arguments and return two values, when the core procedure is
@@ -4672,6 +4712,23 @@ SCM_DEFINE (scm_integer_expt, "integer-expt", 2, 0, 0,
       else  /* return NaN for (0 ^ k) for negative k per R6RS */
        return scm_nan ();
     }
+  /* Optimize the fraction case, to avoid needless reduction of
+     intermediate products to lowest terms.  There will never be any
+     common factor to cancel out, so it would be a waste.  */
+  else if (SCM_FRACTIONP (n))
+    {
+      if (scm_is_true (scm_positive_p (k)))
+       return scm_i_make_ratio_already_reduced
+         (scm_integer_expt (SCM_FRACTION_NUMERATOR (n), k),
+          scm_integer_expt (SCM_FRACTION_DENOMINATOR (n), k));
+      else
+       {
+         k = scm_difference (k, SCM_UNDEFINED);
+         return scm_i_make_ratio_already_reduced
+           (scm_integer_expt (SCM_FRACTION_DENOMINATOR (n), k),
+            scm_integer_expt (SCM_FRACTION_NUMERATOR (n), k));
+       }
+    }
 
   if (SCM_I_INUMP (k))
     i2 = SCM_I_INUM (k);
@@ -7327,8 +7384,9 @@ scm_difference (SCM x, SCM y)
           return scm_c_make_rectangular (-SCM_COMPLEX_REAL (x),
                                    -SCM_COMPLEX_IMAG (x));
        else if (SCM_FRACTIONP (x))
-         return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (x), 
SCM_UNDEFINED),
-                                SCM_FRACTION_DENOMINATOR (x));
+         return scm_i_make_ratio_already_reduced
+           (scm_difference (SCM_FRACTION_NUMERATOR (x), SCM_UNDEFINED),
+            SCM_FRACTION_DENOMINATOR (x));
         else
           SCM_WTA_DISPATCH_1 (g_difference, x, SCM_ARG1, s_difference);
     }
@@ -7876,14 +7934,14 @@ do_divide (SCM x, SCM y, int inexact)
            {
              if (inexact)
                return scm_from_double (1.0 / (double) xx);
-             else return scm_i_make_ratio (SCM_INUM1, x);
+             else return scm_i_make_ratio_already_reduced (SCM_INUM1, x);
            }
        }
       else if (SCM_BIGP (x))
        {
          if (inexact)
            return scm_from_double (1.0 / scm_i_big2dbl (x));
-         else return scm_i_make_ratio (SCM_INUM1, x);
+         else return scm_i_make_ratio_already_reduced (SCM_INUM1, x);
        }
       else if (SCM_REALP (x))
        {
@@ -7913,8 +7971,8 @@ do_divide (SCM x, SCM y, int inexact)
            }
        }
       else if (SCM_FRACTIONP (x))
-       return scm_i_make_ratio (SCM_FRACTION_DENOMINATOR (x),
-                              SCM_FRACTION_NUMERATOR (x));
+       return scm_i_make_ratio_already_reduced (SCM_FRACTION_DENOMINATOR (x),
+                                                SCM_FRACTION_NUMERATOR (x));
       else
        SCM_WTA_DISPATCH_1 (g_divide, x, SCM_ARG1, s_divide);
     }
@@ -8877,8 +8935,9 @@ SCM_PRIMITIVE_GENERIC (scm_magnitude, "magnitude", 1, 0, 
0,
     {
       if (scm_is_false (scm_negative_p (SCM_FRACTION_NUMERATOR (z))))
        return z;
-      return scm_i_make_ratio (scm_difference (SCM_FRACTION_NUMERATOR (z), 
SCM_UNDEFINED),
-                            SCM_FRACTION_DENOMINATOR (z));
+      return scm_i_make_ratio_already_reduced
+       (scm_difference (SCM_FRACTION_NUMERATOR (z), SCM_UNDEFINED),
+        SCM_FRACTION_DENOMINATOR (z));
     }
   else
     SCM_WTA_DISPATCH_1 (g_scm_magnitude, z, SCM_ARG1, s_scm_magnitude);
@@ -8979,8 +9038,9 @@ SCM_PRIMITIVE_GENERIC (scm_inexact_to_exact, 
"inexact->exact", 1, 0, 0,
          
          mpq_init (frac);
          mpq_set_d (frac, val);
-         q = scm_i_make_ratio (scm_i_mpz2num (mpq_numref (frac)),
-                               scm_i_mpz2num (mpq_denref (frac)));
+         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...
diff --git a/test-suite/tests/numbers.test b/test-suite/tests/numbers.test
index 66aa01a..20b7eb1 100644
--- a/test-suite/tests/numbers.test
+++ b/test-suite/tests/numbers.test
@@ -4044,6 +4044,9 @@
   (pass-if (eqv? -0.125 (expt -2 -3.0)))
   (pass-if (eqv? -0.125 (expt -2.0 -3.0)))
   (pass-if (eqv? 0.25 (expt 2.0 -2.0)))
+  (pass-if (eqv? 32/243 (expt 2/3 5)))
+  (pass-if (eqv? 243/32 (expt 2/3 -5)))
+  (pass-if (eqv? 32 (expt 1/2 -5)))
   (pass-if (test-eqv? (* -1.0+0.0i 12398 12398) (expt +12398i 2.0)))
   (pass-if (eqv-loosely? +i (expt -1 0.5)))
   (pass-if (eqv-loosely? +i (expt -1 1/2)))
@@ -4317,6 +4320,9 @@
   (pass-if (eqv? -1/8 (integer-expt -2 -3)))
   (pass-if (eqv? -0.125 (integer-expt -2.0 -3)))
   (pass-if (eqv? 0.25 (integer-expt 2.0 -2)))
+  (pass-if (eqv? 32/243 (integer-expt 2/3 5)))
+  (pass-if (eqv? 243/32 (integer-expt 2/3 -5)))
+  (pass-if (eqv? 32 (integer-expt 1/2 -5)))
   (pass-if (test-eqv? (* -1.0+0.0i 12398 12398) (integer-expt +12398.0i 2))))
 
 
-- 
1.7.10.4

>From 09b12cb11a64a3cecc15fd54a8e3bb20ead0845c Mon Sep 17 00:00:00 2001
From: Mark H Weaver <address@hidden>
Date: Sun, 3 Mar 2013 04:35:09 -0500
Subject: [PATCH 03/10] Add 'round-ash', a rounding arithmetic shift operator

* libguile/numbers.c (left_shift_exact_integer,
  floor_right_shift_exact_integer, round_right_shift_exact_integer):
  New internal static functions to efficiently compute n * 2^count,
  floor(n / 2^count), and round(n / 2^count), respectively, for
  exact integer N and positive exact integer COUNT.  Used to combine
  the implementations of 'ash' and 'round-ash' with minimal code
  duplication.

  (scm_ash): Reimplement in terms of 'left_shift_exact_integer' and
  'floor_right_shift_exact_integer'.  Note that this function
  efficiently computes floor(n * 2^count).

  (scm_round_ash): New procedure to efficiently compute
  round(n * 2^count), for exact integers N and COUNT.  Implemented
  in terms of 'left_shift_exact_integer' and
  'round_right_shift_exact_integer'.

* libguile/numbers.h: Add prototype for scm_round_ash.  Change the
  name of scm_ash's second formal parameter from 'cnt' to 'count'.

* test-suite/tests/numbers.test (round-ash, ash): Add new unified
  testing framework for 'ash' and 'round-ash'.  Previously, the tests
  for 'ash' were not very comprehensive; for example, they did not
  include a single test where the number to be shifted was a bignum.

* doc/ref/api-data.texi (Bitwise Operations): Add documentation for
  'round-ash'.  Improve documentation for `ash'.
---
 doc/ref/api-data.texi         |   44 +++++---
 libguile/numbers.c            |  228 ++++++++++++++++++++++++++++-------------
 libguile/numbers.h            |    3 +-
 test-suite/tests/numbers.test |  114 +++++++++------------
 4 files changed, 237 insertions(+), 152 deletions(-)

diff --git a/doc/ref/api-data.texi b/doc/ref/api-data.texi
index 9da17d8..b33270c 100644
--- a/doc/ref/api-data.texi
+++ b/doc/ref/api-data.texi
@@ -1686,19 +1686,16 @@ starts from 0 for the least significant bit.
 @end lisp
 @end deffn
 
address@hidden {Scheme Procedure} ash n cnt
address@hidden {C Function} scm_ash (n, cnt)
-Return @var{n} shifted left by @var{cnt} bits, or shifted right if
address@hidden is negative.  This is an ``arithmetic'' shift.
-
-This is effectively a multiplication by @m{2^{cnt}, address@hidden, and
-when @var{cnt} is negative it's a division, rounded towards negative
-infinity.  (Note that this is not the same rounding as @code{quotient}
-does.)
address@hidden {Scheme Procedure} ash n count
address@hidden {C Function} scm_ash (n, count)
+Return @math{floor(@var{n} * address@hidden)}, but computed
+more efficiently.  @var{n} and @var{count} must be exact
+integers.
 
-With @var{n} viewed as an infinite precision twos complement,
address@hidden means a left shift introducing zero bits, or a right shift
-dropping bits.
+With @var{n} viewed as an infinite-precision twos-complement
+integer, @code{ash} means a left shift introducing zero bits
+when @var{count} is positive, or a right shift dropping bits
+when @var{count} is negative.  This is an ``arithmetic'' shift.
 
 @lisp
 (number->string (ash #b1 3) 2)     @result{} "1000"
@@ -1709,6 +1706,29 @@ dropping bits.
 @end lisp
 @end deffn
 
address@hidden {Scheme Procedure} round-ash n count
address@hidden {C Function} scm_round_ash (n, count)
+Return @math{round(@var{n} * address@hidden)}, but computed
+more efficiently.  @var{n} and @var{count} must be exact
+integers.
+
+With @var{n} viewed as an infinite-precision twos-complement
+integer, @code{round-ash} means a left shift introducing zero
+bits when @var{count} is positive, or a right shift rounding
+to the nearest integer (with ties going to the nearest even
+integer) when @var{count} is negative.  This is a rounded
+``arithmetic'' shift.
+
address@hidden
+(number->string (round-ash #b1 3) 2)     @result{} \"1000\"
+(number->string (round-ash #b1010 -1) 2) @result{} \"101\"
+(number->string (round-ash #b1010 -2) 2) @result{} \"10\"
+(number->string (round-ash #b1011 -2) 2) @result{} \"11\"
+(number->string (round-ash #b1101 -2) 2) @result{} \"11\"
+(number->string (round-ash #b1110 -2) 2) @result{} \"100\"
address@hidden lisp
address@hidden deffn
+
 @deffn {Scheme Procedure} logcount n
 @deffnx {C Function} scm_logcount (n)
 Return the number of bits in integer @var{n}.  If @var{n} is
diff --git a/libguile/numbers.c b/libguile/numbers.c
index b5ba1a0..aa63fbd 100644
--- a/libguile/numbers.c
+++ b/libguile/numbers.c
@@ -4786,19 +4786,120 @@ SCM_DEFINE (scm_integer_expt, "integer-expt", 2, 0, 0,
 }
 #undef FUNC_NAME
 
+/* n must be an exact integer, and count > 0.
+   Returns n * 2^count. */
+static SCM
+left_shift_exact_integer (SCM n, long count)
+{
+  if (SCM_I_INUMP (n))
+    {
+      scm_t_inum nn = SCM_I_INUM (n);
+
+      /* Left shift of count >= SCM_I_FIXNUM_BIT-1 will always
+         overflow a non-zero fixnum.  For smaller shifts we check the
+         bits going into positions above SCM_I_FIXNUM_BIT-1.  If they're
+         all 0s for nn>=0, or all 1s for nn<0 then there's no overflow.
+         Those bits are "nn >> (SCM_I_FIXNUM_BIT-1 - count)".  */
+
+      if (nn == 0)
+        return n;
+      else if (count < SCM_I_FIXNUM_BIT-1 &&
+               ((scm_t_bits) (SCM_SRS (nn, (SCM_I_FIXNUM_BIT-1 - count)) + 1)
+                <= 1))
+        return SCM_I_MAKINUM (nn << count);
+      else
+        {
+          SCM result = scm_i_inum2big (nn);
+          mpz_mul_2exp (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (result),
+                        count);
+          return result;
+        }
+    }
+  else if (SCM_BIGP (n))
+    {
+      SCM result = scm_i_mkbig ();
+      mpz_mul_2exp (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (n), count);
+      scm_remember_upto_here_1 (n);
+      return result;
+    }
+  else
+    scm_syserror ("left_shift_exact_integer");
+}
+
+/* n must be an exact integer, and count > 0.
+   Returns floor(n / 2^count). */
+static SCM
+floor_right_shift_exact_integer (SCM n, long count)
+{
+  if (SCM_I_INUMP (n))
+    {
+      scm_t_inum nn = SCM_I_INUM (n);
+
+      if (count >= SCM_I_FIXNUM_BIT)
+        return (nn >= 0 ? SCM_INUM0 : SCM_I_MAKINUM (-1));
+      else
+        return SCM_I_MAKINUM (SCM_SRS (nn, count));
+    }
+  else if (SCM_BIGP (n))
+    {
+      SCM result = scm_i_mkbig ();
+      mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (n),
+                       count);
+      scm_remember_upto_here_1 (n);
+      return scm_i_normbig (result);
+    }
+  else
+    scm_syserror ("floor_right_shift_exact_integer");
+}
+
+/* n must be an exact integer, and count > 0.
+   Returns round(n / 2^count). */
+static SCM
+round_right_shift_exact_integer (SCM n, long count)
+{
+  if (SCM_I_INUMP (n))
+    {
+      if (count >= SCM_I_FIXNUM_BIT)
+        return SCM_INUM0;
+      else
+        {
+          scm_t_inum nn = SCM_I_INUM (n);
+          scm_t_inum qq = SCM_SRS (nn, count);
+
+          if (0 == (nn & (1L << (count-1))))
+            return SCM_I_MAKINUM (qq);                /* round down */
+          else if (nn & ((1L << (count-1)) - 1))
+            return SCM_I_MAKINUM (qq + 1);            /* round up */
+          else
+            return SCM_I_MAKINUM ((~1L) & (qq + 1));  /* round to even */
+        }
+    }
+  else if (SCM_BIGP (n))
+    {
+      SCM q = scm_i_mkbig ();
+
+      mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (n), count);
+      if (mpz_tstbit (SCM_I_BIG_MPZ (n), count-1)
+          && (mpz_odd_p (SCM_I_BIG_MPZ (q))
+              || (mpz_scan1 (SCM_I_BIG_MPZ (n), 0) < count-1)))
+        mpz_add_ui (SCM_I_BIG_MPZ (q), SCM_I_BIG_MPZ (q), 1);
+      scm_remember_upto_here_1 (n);
+      return scm_i_normbig (q);
+    }
+  else
+    scm_syserror ("round_right_shift_exact_integer");
+}
+
 SCM_DEFINE (scm_ash, "ash", 2, 0, 0,
-            (SCM n, SCM cnt),
-           "Return @var{n} shifted left by @var{cnt} bits, or shifted right\n"
-           "if @var{cnt} is negative.  This is an ``arithmetic'' shift.\n"
+            (SCM n, SCM count),
+           "Return @math{floor(@var{n} * address@hidden)}, but computed\n"
+           "more efficiently.  @var{n} and @var{count} must be exact\n"
+           "integers.\n"
            "\n"
-           "This is effectively a multiplication by address@hidden, and when\n"
-           "@var{cnt} is negative it's a division, rounded towards negative\n"
-           "infinity.  (Note that this is not the same rounding as\n"
-           "@code{quotient} does.)\n"
-           "\n"
-           "With @var{n} viewed as an infinite precision twos complement,\n"
-           "@code{ash} means a left shift introducing zero bits, or a right\n"
-           "shift dropping bits.\n"
+           "With @var{n} viewed as an infinite-precision twos-complement\n"
+           "integer, @code{ash} means a left shift introducing zero bits\n"
+           "when @var{count} is positive, or a right shift dropping bits\n"
+           "when @var{count} is negative.  This is an ``arithmetic'' shift.\n"
            "\n"
            "@lisp\n"
            "(number->string (ash #b1 3) 2)     @result{} \"1000\"\n"
@@ -4809,79 +4910,58 @@ SCM_DEFINE (scm_ash, "ash", 2, 0, 0,
            "@end lisp")
 #define FUNC_NAME s_scm_ash
 {
-  long bits_to_shift;
-  bits_to_shift = scm_to_long (cnt);
-
-  if (SCM_I_INUMP (n))
+  if (SCM_I_INUMP (n) || SCM_BIGP (n))
     {
-      scm_t_inum nn = SCM_I_INUM (n);
+      long bits_to_shift = scm_to_long (count);
 
       if (bits_to_shift > 0)
-        {
-          /* Left shift of bits_to_shift >= SCM_I_FIXNUM_BIT-1 will always
-             overflow a non-zero fixnum.  For smaller shifts we check the
-             bits going into positions above SCM_I_FIXNUM_BIT-1.  If they're
-             all 0s for nn>=0, or all 1s for nn<0 then there's no overflow.
-             Those bits are "nn >> (SCM_I_FIXNUM_BIT-1 -
-             bits_to_shift)".  */
-
-          if (nn == 0)
-            return n;
-
-          if (bits_to_shift < SCM_I_FIXNUM_BIT-1
-              && ((scm_t_bits)
-                  (SCM_SRS (nn, (SCM_I_FIXNUM_BIT-1 - bits_to_shift)) + 1)
-                  <= 1))
-            {
-              return SCM_I_MAKINUM (nn << bits_to_shift);
-            }
-          else
-            {
-              SCM result = scm_i_inum2big (nn);
-              mpz_mul_2exp (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (result),
-                            bits_to_shift);
-              return result;
-            }
-        }
+        return left_shift_exact_integer (n, bits_to_shift);
+      else if (SCM_LIKELY (bits_to_shift < 0))
+        return floor_right_shift_exact_integer (n, -bits_to_shift);
       else
-        {
-          bits_to_shift = -bits_to_shift;
-          if (bits_to_shift >= SCM_LONG_BIT)
-            return (nn >= 0 ? SCM_INUM0 : SCM_I_MAKINUM(-1));
-          else
-            return SCM_I_MAKINUM (SCM_SRS (nn, bits_to_shift));
-        }
-
+        return n;
     }
-  else if (SCM_BIGP (n))
-    {
-      SCM result;
+  else
+    SCM_WRONG_TYPE_ARG (SCM_ARG1, n);
+}
+#undef FUNC_NAME
 
-      if (bits_to_shift == 0)
-        return n;
+SCM_DEFINE (scm_round_ash, "round-ash", 2, 0, 0,
+            (SCM n, SCM count),
+           "Return @math{round(@var{n} * address@hidden)}, but computed\n"
+           "more efficiently.  @var{n} and @var{count} must be exact\n"
+           "integers.\n"
+           "\n"
+           "With @var{n} viewed as an infinite-precision twos-complement\n"
+           "integer, @code{round-ash} means a left shift introducing zero\n"
+           "bits when @var{count} is positive, or a right shift rounding\n"
+           "to the nearest integer (with ties going to the nearest even\n"
+           "integer) when @var{count} is negative.  This is a rounded\n"
+           "``arithmetic'' shift.\n"
+           "\n"
+           "@lisp\n"
+           "(number->string (round-ash #b1 3) 2)     @result{} \"1000\"\n"
+           "(number->string (round-ash #b1010 -1) 2) @result{} \"101\"\n"
+           "(number->string (round-ash #b1010 -2) 2) @result{} \"10\"\n"
+           "(number->string (round-ash #b1011 -2) 2) @result{} \"11\"\n"
+           "(number->string (round-ash #b1101 -2) 2) @result{} \"11\"\n"
+           "(number->string (round-ash #b1110 -2) 2) @result{} \"100\"\n"
+           "@end lisp")
+#define FUNC_NAME s_scm_round_ash
+{
+  if (SCM_I_INUMP (n) || SCM_BIGP (n))
+    {
+      long bits_to_shift = scm_to_long (count);
 
-      result = scm_i_mkbig ();
-      if (bits_to_shift >= 0)
-        {
-          mpz_mul_2exp (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (n),
-                        bits_to_shift);
-          return result;
-        }
+      if (bits_to_shift > 0)
+        return left_shift_exact_integer (n, bits_to_shift);
+      else if (SCM_LIKELY (bits_to_shift < 0))
+        return round_right_shift_exact_integer (n, -bits_to_shift);
       else
-        {
-          /* GMP doesn't have an fdiv_q_2exp variant returning just a long, so
-             we have to allocate a bignum even if the result is going to be a
-             fixnum.  */
-          mpz_fdiv_q_2exp (SCM_I_BIG_MPZ (result), SCM_I_BIG_MPZ (n),
-                           -bits_to_shift);
-          return scm_i_normbig (result);
-        }
-
+        return n;
     }
   else
-    {
-      SCM_WRONG_TYPE_ARG (SCM_ARG1, n);
-    }
+    SCM_WRONG_TYPE_ARG (SCM_ARG1, n);
 }
 #undef FUNC_NAME
 
diff --git a/libguile/numbers.h b/libguile/numbers.h
index 2c8b260..912f287 100644
--- a/libguile/numbers.h
+++ b/libguile/numbers.h
@@ -206,7 +206,8 @@ SCM_API SCM scm_logbit_p (SCM n1, SCM n2);
 SCM_API SCM scm_lognot (SCM n);
 SCM_API SCM scm_modulo_expt (SCM n, SCM k, SCM m);
 SCM_API SCM scm_integer_expt (SCM z1, SCM z2);
-SCM_API SCM scm_ash (SCM n, SCM cnt);
+SCM_API SCM scm_ash (SCM n, SCM count);
+SCM_API SCM scm_round_ash (SCM n, SCM count);
 SCM_API SCM scm_bit_extract (SCM n, SCM start, SCM end);
 SCM_API SCM scm_logcount (SCM n);
 SCM_API SCM scm_integer_length (SCM n);
diff --git a/test-suite/tests/numbers.test b/test-suite/tests/numbers.test
index 20b7eb1..8dab29c 100644
--- a/test-suite/tests/numbers.test
+++ b/test-suite/tests/numbers.test
@@ -201,71 +201,6 @@
     (eqv? -2305843009213693953 (1- -2305843009213693952))))
 
 ;;;
-;;; ash
-;;;
-
-(with-test-prefix "ash"
-
-  (pass-if "documented?"
-    (documented? ash))
-
-  (pass-if (eqv? 0 (ash 0 0)))
-  (pass-if (eqv? 0 (ash 0 1)))
-  (pass-if (eqv? 0 (ash 0 1000)))
-  (pass-if (eqv? 0 (ash 0 -1)))
-  (pass-if (eqv? 0 (ash 0 -1000)))
-
-  (pass-if (eqv? 1 (ash 1 0)))
-  (pass-if (eqv? 2 (ash 1 1)))
-  (pass-if (eqv? 340282366920938463463374607431768211456 (ash 1 128)))
-  (pass-if (eqv? 0 (ash 1 -1)))
-  (pass-if (eqv? 0 (ash 1 -1000)))
-
-  (pass-if (eqv? -1 (ash -1 0)))
-  (pass-if (eqv? -2 (ash -1 1)))
-  (pass-if (eqv? -340282366920938463463374607431768211456 (ash -1 128)))
-  (pass-if (eqv? -1 (ash -1 -1)))
-  (pass-if (eqv? -1 (ash -1 -1000)))
-
-  (pass-if (eqv? -3 (ash -3 0)))
-  (pass-if (eqv? -6 (ash -3 1)))
-  (pass-if (eqv? -1020847100762815390390123822295304634368 (ash -3 128)))
-  (pass-if (eqv? -2 (ash -3 -1)))
-  (pass-if (eqv? -1 (ash -3 -1000)))
-
-  (pass-if (eqv? -6 (ash -23 -2)))
-
-  (pass-if (eqv? most-positive-fixnum       (ash most-positive-fixnum 0)))
-  (pass-if (eqv? (* 2 most-positive-fixnum) (ash most-positive-fixnum 1)))
-  (pass-if (eqv? (* 4 most-positive-fixnum) (ash most-positive-fixnum 2)))
-  (pass-if
-      (eqv? (* most-positive-fixnum 340282366920938463463374607431768211456)
-           (ash most-positive-fixnum 128)))
-  (pass-if (eqv? (quotient most-positive-fixnum 2)
-                (ash most-positive-fixnum -1)))
-  (pass-if (eqv? 0 (ash most-positive-fixnum -1000)))
-
-  (let ((mpf4 (quotient most-positive-fixnum 4)))
-    (pass-if (eqv? (* 2 mpf4) (ash mpf4 1)))
-    (pass-if (eqv? (* 4 mpf4) (ash mpf4 2)))
-    (pass-if (eqv? (* 8 mpf4) (ash mpf4 3))))
-
-  (pass-if (eqv? most-negative-fixnum       (ash most-negative-fixnum 0)))
-  (pass-if (eqv? (* 2 most-negative-fixnum) (ash most-negative-fixnum 1)))
-  (pass-if (eqv? (* 4 most-negative-fixnum) (ash most-negative-fixnum 2)))
-  (pass-if
-      (eqv? (* most-negative-fixnum 340282366920938463463374607431768211456)
-           (ash most-negative-fixnum 128)))
-  (pass-if (eqv? (quotient-floor most-negative-fixnum 2)
-                (ash most-negative-fixnum -1)))
-  (pass-if (eqv? -1 (ash most-negative-fixnum -1000)))
-
-  (let ((mnf4 (quotient-floor most-negative-fixnum 4)))
-    (pass-if (eqv? (* 2 mnf4) (ash mnf4 1)))
-    (pass-if (eqv? (* 4 mnf4) (ash mnf4 2)))
-    (pass-if (eqv? (* 8 mnf4) (ash mnf4 3)))))
-
-;;;
 ;;; exact?
 ;;;
 
@@ -4904,3 +4839,52 @@
                         round-quotient
                         round-remainder
                         valid-round-answer?)))
+
+;;;
+;;; ash
+;;; round-ash
+;;;
+
+(let ()
+  (define (test-ash-variant name ash-variant round-variant)
+    (with-test-prefix name
+      (define (test n count)
+        (pass-if (list n count)
+          (eqv? (ash-variant n count)
+                (round-variant (* n (expt 2 count))))))
+
+      (pass-if "documented?"
+        (documented? ash-variant))
+
+      (for-each (lambda (n)
+                  (for-each (lambda (count) (test n count))
+                            '(-1000 -3 -2 -1 0 1 2 3 1000)))
+                (list 0 1 3 23 -1 -3 -23
+                      fixnum-max
+                      (1+ fixnum-max)
+                      (1- fixnum-max)
+                      (* fixnum-max 4)
+                      (quotient fixnum-max 4)
+                      fixnum-min
+                      (1+ fixnum-min)
+                      (1- fixnum-min)
+                      (* fixnum-min 4)
+                      (quotient fixnum-min 4)))
+
+      (do ((count -2 (1- count))
+           (vals '(1 3 5 7 9 11)
+                 (map (lambda (n) (* 2 n)) vals)))
+          ((> (car vals) (* 2 fixnum-max)) 'done)
+        (for-each (lambda (n)
+                    (test    n  count)
+                    (test (- n) count))
+                  vals))
+
+      ;; Test rounding
+      (for-each (lambda (base)
+                  (for-each (lambda (offset) (test (+ base offset) -3))
+                            '(#b11001 #b11100 #b11101 #b10001 #b10100 
#b10101)))
+                (list 0 64 -64 (* 64 fixnum-max) (* 64 fixnum-min)))))
+
+  (test-ash-variant       'ash       ash floor)
+  (test-ash-variant 'round-ash round-ash round))
-- 
1.7.10.4

>From e163f7fbf1b3af95fef5ca13f8a36166cbdb0233 Mon Sep 17 00:00:00 2001
From: Mark H Weaver <address@hidden>
Date: Mon, 4 Mar 2013 18:37:23 -0500
Subject: [PATCH 04/10] Verify that FLT_RADIX is 2.

* libguile/numbers.c: Trigger a compilation error if FLT_RADIX is not 2.
  This has long been assumed by code in numbers.c.
---
 libguile/numbers.c |    3 +++
 1 file changed, 3 insertions(+)

diff --git a/libguile/numbers.c b/libguile/numbers.c
index aa63fbd..f9c65a8 100644
--- a/libguile/numbers.c
+++ b/libguile/numbers.c
@@ -81,6 +81,9 @@
 #define M_PI       3.14159265358979323846
 #endif
 
+/* FIXME: We assume that FLT_RADIX is 2 */
+verify (FLT_RADIX == 2);
+
 typedef scm_t_signed_bits scm_t_inum;
 #define scm_from_inum(x) (scm_from_signed_integer (x))
 
-- 
1.7.10.4

>From ccd155332040dff20dff0dde63fd2a7656f0900a Mon Sep 17 00:00:00 2001
From: Mark H Weaver <address@hidden>
Date: Sun, 3 Mar 2013 04:58:55 -0500
Subject: [PATCH 05/10] Simplify and improve scm_i_big2dbl, and add
 scm_i_big2dbl_2exp

* libguile/numbers.c: Verify that FLT_RADIX == 2, which has long
  been assumed in various places.

  (scm_i_big2dbl_2exp): New static function.

  (scm_i_big2dbl): Reimplement in terms of 'scm_i_big2dbl_2exp'.
  Ties now go to even.
---
 libguile/numbers.c |  101 +++++++++++++++++++---------------------------------
 1 file changed, 36 insertions(+), 65 deletions(-)

diff --git a/libguile/numbers.c b/libguile/numbers.c
index f9c65a8..11911d7 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
-- 
1.7.10.4

>From 5fc88dd6f30a693eff5bea2faadb69c12d3a8359 Mon Sep 17 00:00:00 2001
From: Mark H Weaver <address@hidden>
Date: Sun, 3 Mar 2013 05:02:53 -0500
Subject: [PATCH 06/10] Optimize logarithms using scm_i_big2dbl_2exp

* libguile/numbers.c (log_of_exact_integer_with_size): Removed.

  (log_of_exact_integer): Handle bignums too large to fit in a double
  using scm_i_big2dbl_2exp instead of scm_integer_length and scm_ash.

  (log_of_fraction): Use log_of_exact_integer instead of
  log_of_exact_integer_with_size.
---
 libguile/numbers.c |   30 ++++++++++++------------------
 1 file changed, 12 insertions(+), 18 deletions(-)

diff --git a/libguile/numbers.c b/libguile/numbers.c
index 11911d7..e9059a3 100644
--- a/libguile/numbers.c
+++ b/libguile/numbers.c
@@ -9576,26 +9576,20 @@ log_of_shifted_double (double x, long shift)
     return scm_c_make_rectangular (ans, M_PI);
 }
 
-/* Returns log(n), for exact integer n of integer-length size */
-static SCM
-log_of_exact_integer_with_size (SCM n, long size)
-{
-  long shift = size - 2 * scm_dblprec[0];
-
-  if (shift > 0)
-    return log_of_shifted_double
-      (scm_to_double (scm_ash (n, scm_from_long(-shift))),
-       shift);
-  else
-    return log_of_shifted_double (scm_to_double (n), 0);
-}
-
 /* Returns log(n), for exact integer n */
 static SCM
 log_of_exact_integer (SCM n)
 {
-  return log_of_exact_integer_with_size
-    (n, scm_to_long (scm_integer_length (n)));
+  if (SCM_I_INUMP (n))
+    return log_of_shifted_double (SCM_I_INUM (n), 0);
+  else if (SCM_BIGP (n))
+    {
+      long expon;
+      double signif = scm_i_big2dbl_2exp (n, &expon);
+      return log_of_shifted_double (signif, expon);
+    }
+  else
+    scm_wrong_type_arg ("log_of_exact_integer", SCM_ARG1, n);
 }
 
 /* Returns log(n/d), for exact non-zero integers n and d */
@@ -9606,8 +9600,8 @@ log_of_fraction (SCM n, SCM d)
   long d_size = scm_to_long (scm_integer_length (d));
 
   if (abs (n_size - d_size) > 1)
-    return (scm_difference (log_of_exact_integer_with_size (n, n_size),
-                           log_of_exact_integer_with_size (d, d_size)));
+    return (scm_difference (log_of_exact_integer (n),
+                           log_of_exact_integer (d)));
   else if (scm_is_false (scm_negative_p (n)))
     return scm_from_double
       (log1p (scm_to_double (scm_divide2real (scm_difference (n, d), d))));
-- 
1.7.10.4

>From fcd504b656a47ef0b92d916c2503214ed0220002 Mon Sep 17 00:00:00 2001
From: Mark H Weaver <address@hidden>
Date: Mon, 4 Mar 2013 18:42:27 -0500
Subject: [PATCH 07/10] 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.
---
 libguile/numbers.c |   41 +++++++++++++++++++++++++++--------------
 1 file changed, 27 insertions(+), 14 deletions(-)

diff --git a/libguile/numbers.c b/libguile/numbers.c
index e9059a3..5feed70 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;
        }
     }
 }
-- 
1.7.10.4

>From 3f5f479b73744d939d15aa1a5736e013ded3ea1d Mon Sep 17 00:00:00 2001
From: Mark H Weaver <address@hidden>
Date: Mon, 4 Mar 2013 18:46:33 -0500
Subject: [PATCH 08/10] Improve inexact division of exact integers.

* libguile/numbers.c (scm_i_divide2double): New function.
  (scm_i_divide2double_lo2b): New variable.
  (scm_i_fraction2double, log_of_fraction): Use 'scm_i_divide2double'.
  (do_divide): Removed.  Its code is now in 'scm_divide'.
  (scm_divide2real): Removed.  Superceded by 'scm_i_divide2double'.
  (scm_divide): Inherit code from 'do_divide', but without support for
  forcing a 'double' result (that functionality is now implemented by
  'scm_i_divide2double').  Add FIXME comments in cases where divisions
  might not be as precise as they should be.
  (scm_init_numbers): Initialize 'scm_i_divide2double_lo2b'.
---
 libguile/numbers.c |  257 ++++++++++++++++++++++++++++++++++++----------------
 1 file changed, 177 insertions(+), 80 deletions(-)

diff --git a/libguile/numbers.c b/libguile/numbers.c
index 5feed70..601b2a6 100644
--- a/libguile/numbers.c
+++ b/libguile/numbers.c
@@ -410,9 +410,6 @@ scm_i_mpz2num (mpz_t b)
   }
 }
 
-/* this is needed when we want scm_divide to make a float, not a ratio, even 
if passed two ints */
-static SCM scm_divide2real (SCM x, SCM y);
-
 /* scm_i_make_ratio_already_reduced makes a ratio more efficiently,
    when the following conditions are known in advance:
 
@@ -468,11 +465,122 @@ scm_i_make_ratio (SCM numerator, SCM denominator)
 }
 #undef FUNC_NAME
 
+static mpz_t scm_i_divide2double_lo2b;
+
+/* Return n/d as a double, where n and d are exact integers.  */
+static double
+scm_i_divide2double (SCM n, SCM d)
+{
+  int neg;
+  mpz_t nn, dd, lo, hi, x;
+  ssize_t e;
+
+  if (SCM_I_INUMP (d))
+    {
+      if (SCM_UNLIKELY (scm_is_eq (d, SCM_INUM0)))
+        {
+          if (scm_is_true (scm_positive_p (n)))
+            return 1.0 / 0.0;
+          else if (scm_is_true (scm_negative_p (n)))
+            return -1.0 / 0.0;
+          else
+            return 0.0 / 0.0;
+        }
+      mpz_init_set_si (dd, SCM_I_INUM (d));
+    }
+  else
+    mpz_init_set (dd, SCM_I_BIG_MPZ (d));
+
+  if (SCM_I_INUMP (n))
+    mpz_init_set_si (nn, SCM_I_INUM (n));
+  else
+    mpz_init_set (nn, SCM_I_BIG_MPZ (n));
+
+  neg = (mpz_sgn (nn) < 0) ^ (mpz_sgn (dd) < 0);
+  mpz_abs (nn, nn);
+  mpz_abs (dd, dd);
+
+  /* Now we need to find the value of e such that:
+
+     For e <= 0:
+       b^{p-1} - 1/2b <= b^-e n / d < b^p - 1/2
+       (2 b^p - 1) <= 2 b b^-e n / d < (2 b^p - 1) b
+       (2 b^p - 1) d <= 2 b b^-e n < (2 b^p - 1) d b
+
+     For e >= 0:
+       b^{p-1} - 1/2b <= n / b^e d < b^p - 1/2    
+       (2 b^p - 1) <= 2 b n / b^e d < (2 b^p - 1) b
+       (2 b^p - 1) d b^e <= 2 b n < (2 b^p - 1) d b b^e
+
+     p == DBL_MANT_DIG
+     b == FLT_RADIX (assumed 2 here) */
+
+  /* Initial guess for e */
+  e = mpz_sizeinbase (nn, 2) - mpz_sizeinbase (dd, 2) - (DBL_MANT_DIG-1);
+  if (e < DBL_MIN_EXP - DBL_MANT_DIG)
+    e = DBL_MIN_EXP - DBL_MANT_DIG;
+
+  /* Compute initial lo, hi, and x */
+  mpz_inits (lo, hi, x, NULL);
+  mpz_mul_2exp (x, nn, 2 + ((e < 0) ? -e : 0));
+  mpz_mul (lo, dd, scm_i_divide2double_lo2b);
+  if (e > 0)
+    mpz_mul_2exp (lo, lo, e);
+  mpz_mul_2exp (hi, lo, 1);
+
+  /* Adjust e as needed */
+  while (mpz_cmp (x, lo) < 0 && e > DBL_MIN_EXP - DBL_MANT_DIG)
+    {
+      mpz_mul_2exp (x, x, 1);
+      e--;
+    }
+  while (mpz_cmp (x, hi) >= 0)
+    {
+      mpz_mul_2exp (hi, hi, 1);
+      e++;
+    }
+
+  /* Now compute the rounded mantissa:
+     n / b^e d   (if e >= 0)
+     n b^-e / d  (if e <= 0) */
+  {
+    int cmp, needs_adjustment;
+    double result;
+
+    if (e < 0)
+      mpz_mul_2exp (nn, nn, -e);
+    else
+      mpz_mul_2exp (dd, dd, e);
+
+    /* mpz does not directly support rounded right
+       shifts, so we have to do it the hard way.
+       hi == quotient, lo == remainder */
+    mpz_fdiv_qr (hi, lo, nn, dd);
+    mpz_mul_2exp (lo, lo, 1);
+
+    cmp = mpz_cmp (lo, dd);
+    if (mpz_odd_p (hi))
+      needs_adjustment = (cmp >= 0);
+    else
+      needs_adjustment = (cmp > 0);
+
+    if (needs_adjustment)
+      mpz_add_ui (hi, hi, 1);
+
+    if (neg)
+      mpz_neg (hi, hi);
+
+    result = ldexp (mpz_get_d (hi), e);
+    mpz_clears (nn, dd, lo, hi, x, NULL);
+    return result;
+  }
+}
+
 double
 scm_i_fraction2double (SCM z)
 {
-  return scm_to_double (scm_divide2real (SCM_FRACTION_NUMERATOR (z), 
-                                        SCM_FRACTION_DENOMINATOR (z)));
+  return scm_i_divide2double (SCM_FRACTION_NUMERATOR (z),
+                              SCM_FRACTION_DENOMINATOR (z));
 }
 
 static int
@@ -7965,8 +8073,8 @@ SCM_PRIMITIVE_GENERIC (scm_i_divide, "/", 0, 2, 1,
 #define s_divide s_scm_i_divide
 #define g_divide g_scm_i_divide
 
-static SCM
-do_divide (SCM x, SCM y, int inexact)
+SCM
+scm_divide (SCM x, SCM y)
 #define FUNC_NAME s_divide
 {
   double a;
@@ -7985,18 +8093,10 @@ do_divide (SCM x, SCM y, int inexact)
            scm_num_overflow (s_divide);
 #endif
          else
-           {
-             if (inexact)
-               return scm_from_double (1.0 / (double) xx);
-             else return scm_i_make_ratio_already_reduced (SCM_INUM1, x);
-           }
+           return scm_i_make_ratio_already_reduced (SCM_INUM1, x);
        }
       else if (SCM_BIGP (x))
-       {
-         if (inexact)
-           return scm_from_double (1.0 / scm_i_big2dbl (x));
-         else return scm_i_make_ratio_already_reduced (SCM_INUM1, x);
-       }
+       return scm_i_make_ratio_already_reduced (SCM_INUM1, x);
       else if (SCM_REALP (x))
        {
          double xx = SCM_REAL_VALUE (x);
@@ -8046,11 +8146,7 @@ do_divide (SCM x, SCM y, int inexact)
 #endif
            }
          else if (xx % yy != 0)
-           {
-             if (inexact)
-               return scm_from_double ((double) xx / (double) yy);
-             else return scm_i_make_ratio (x, y);
-           }
+           return scm_i_make_ratio (x, y);
          else
            {
              scm_t_inum z = xx / yy;
@@ -8061,11 +8157,7 @@ do_divide (SCM x, SCM y, int inexact)
            }
        }
       else if (SCM_BIGP (y))
-       {
-         if (inexact)
-           return scm_from_double ((double) xx / scm_i_big2dbl (y));
-         else return scm_i_make_ratio (x, y);
-       }
+       return scm_i_make_ratio (x, y);
       else if (SCM_REALP (y))
        {
          double yy = SCM_REAL_VALUE (y);
@@ -8074,6 +8166,9 @@ do_divide (SCM x, SCM y, int inexact)
            scm_num_overflow (s_divide);
          else
 #endif
+            /* FIXME: Precision may be lost here due to:
+               (1) The cast from 'scm_t_inum' to 'double'
+               (2) Double rounding */
            return scm_from_double ((double) xx / yy);
        }
       else if (SCM_COMPLEXP (y))
@@ -8100,7 +8195,7 @@ do_divide (SCM x, SCM y, int inexact)
       else if (SCM_FRACTIONP (y))
        /* a / b/c = ac / b */
        return scm_i_make_ratio (scm_product (x, SCM_FRACTION_DENOMINATOR (y)),
-                              SCM_FRACTION_NUMERATOR (y));
+                                 SCM_FRACTION_NUMERATOR (y));
       else
        SCM_WTA_DISPATCH_2 (g_divide, x, y, SCM_ARGn, s_divide);
     }
@@ -8144,43 +8239,24 @@ do_divide (SCM x, SCM y, int inexact)
                  return scm_i_normbig (result);
                }
              else
-               {
-                 if (inexact)
-                   return scm_from_double (scm_i_big2dbl (x) / (double) yy);
-                 else return scm_i_make_ratio (x, y);
-               }
+               return scm_i_make_ratio (x, y);
            }
        }
       else if (SCM_BIGP (y))
        {
-         /* big_x / big_y */
-         if (inexact)
-           {
-             /* It's easily possible for the ratio x/y to fit a double
-                but one or both x and y be too big to fit a double,
-                hence the use of mpq_get_d rather than converting and
-                dividing.  */
-             mpq_t q;
-             *mpq_numref(q) = *SCM_I_BIG_MPZ (x);
-             *mpq_denref(q) = *SCM_I_BIG_MPZ (y);
-             return scm_from_double (mpq_get_d (q));
-           }
-         else
-           {
-             int divisible_p = mpz_divisible_p (SCM_I_BIG_MPZ (x),
-                                                SCM_I_BIG_MPZ (y));
-             if (divisible_p)
-               {
-                 SCM result = scm_i_mkbig ();
-                 mpz_divexact (SCM_I_BIG_MPZ (result),
-                               SCM_I_BIG_MPZ (x),
-                               SCM_I_BIG_MPZ (y));
-                 scm_remember_upto_here_2 (x, y);
-                 return scm_i_normbig (result);
-               }
-             else
-               return scm_i_make_ratio (x, y);
-           }
+          int divisible_p = mpz_divisible_p (SCM_I_BIG_MPZ (x),
+                                             SCM_I_BIG_MPZ (y));
+          if (divisible_p)
+            {
+              SCM result = scm_i_mkbig ();
+              mpz_divexact (SCM_I_BIG_MPZ (result),
+                            SCM_I_BIG_MPZ (x),
+                            SCM_I_BIG_MPZ (y));
+              scm_remember_upto_here_2 (x, y);
+              return scm_i_normbig (result);
+            }
+          else
+            return scm_i_make_ratio (x, y);
        }
       else if (SCM_REALP (y))
        {
@@ -8190,6 +8266,8 @@ do_divide (SCM x, SCM y, int inexact)
            scm_num_overflow (s_divide);
          else
 #endif
+            /* FIXME: Precision may be lost here due to:
+               (1) scm_i_big2dbl (2) Double rounding */
            return scm_from_double (scm_i_big2dbl (x) / yy);
        }
       else if (SCM_COMPLEXP (y))
@@ -8199,7 +8277,7 @@ do_divide (SCM x, SCM y, int inexact)
        }
       else if (SCM_FRACTIONP (y))
        return scm_i_make_ratio (scm_product (x, SCM_FRACTION_DENOMINATOR (y)),
-                              SCM_FRACTION_NUMERATOR (y));
+                                 SCM_FRACTION_NUMERATOR (y));
       else
        SCM_WTA_DISPATCH_2 (g_divide, x, y, SCM_ARGn, s_divide);
     }
@@ -8214,10 +8292,16 @@ do_divide (SCM x, SCM y, int inexact)
            scm_num_overflow (s_divide);
          else
 #endif
+            /* FIXME: Precision may be lost here due to:
+               (1) The cast from 'scm_t_inum' to 'double'
+               (2) Double rounding */
            return scm_from_double (rx / (double) yy);
        }
       else if (SCM_BIGP (y))
        {
+          /* FIXME: Precision may be lost here due to:
+             (1) The conversion from bignum to double
+             (2) Double rounding */
          double dby = mpz_get_d (SCM_I_BIG_MPZ (y));
          scm_remember_upto_here_1 (y);
          return scm_from_double (rx / dby);
@@ -8255,12 +8339,18 @@ do_divide (SCM x, SCM y, int inexact)
          else
 #endif
            {
+              /* FIXME: Precision may be lost here due to:
+                 (1) The conversion from 'scm_t_inum' to double
+                 (2) Double rounding */
              double d = yy;
              return scm_c_make_rectangular (rx / d, ix / d);
            }
        }
       else if (SCM_BIGP (y))
        {
+          /* FIXME: Precision may be lost here due to:
+             (1) The conversion from bignum to double
+             (2) Double rounding */
          double dby = mpz_get_d (SCM_I_BIG_MPZ (y));
          scm_remember_upto_here_1 (y);
          return scm_c_make_rectangular (rx / dby, ix / dby);
@@ -8294,6 +8384,9 @@ do_divide (SCM x, SCM y, int inexact)
        }
       else if (SCM_FRACTIONP (y))
        {
+          /* FIXME: Precision may be lost here due to:
+             (1) The conversion from fraction to double
+             (2) Double rounding */
          double yy = scm_i_fraction2double (y);
          return scm_c_make_rectangular (rx / yy, ix / yy);
        }
@@ -8311,12 +8404,12 @@ do_divide (SCM x, SCM y, int inexact)
          else
 #endif
            return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x),
-                                  scm_product (SCM_FRACTION_DENOMINATOR (x), 
y));
+                                     scm_product (SCM_FRACTION_DENOMINATOR 
(x), y));
        } 
       else if (SCM_BIGP (y)) 
        {
          return scm_i_make_ratio (SCM_FRACTION_NUMERATOR (x),
-                                scm_product (SCM_FRACTION_DENOMINATOR (x), y));
+                                   scm_product (SCM_FRACTION_DENOMINATOR (x), 
y));
        } 
       else if (SCM_REALP (y)) 
        {
@@ -8326,33 +8419,28 @@ do_divide (SCM x, SCM y, int inexact)
            scm_num_overflow (s_divide);
          else
 #endif
+            /* FIXME: Precision may be lost here due to:
+               (1) The conversion from fraction to double
+               (2) Double rounding */
            return scm_from_double (scm_i_fraction2double (x) / yy);
        }
       else if (SCM_COMPLEXP (y)) 
        {
+          /* FIXME: Precision may be lost here due to:
+             (1) The conversion from fraction to double
+             (2) Double rounding */
          a = scm_i_fraction2double (x);
          goto complex_div;
        } 
       else if (SCM_FRACTIONP (y))
        return scm_i_make_ratio (scm_product (SCM_FRACTION_NUMERATOR (x), 
SCM_FRACTION_DENOMINATOR (y)),
-                              scm_product (SCM_FRACTION_NUMERATOR (y), 
SCM_FRACTION_DENOMINATOR (x)));
+                                 scm_product (SCM_FRACTION_NUMERATOR (y), 
SCM_FRACTION_DENOMINATOR (x)));
       else 
        SCM_WTA_DISPATCH_2 (g_divide, x, y, SCM_ARGn, s_divide);
     }
   else
     SCM_WTA_DISPATCH_2 (g_divide, x, y, SCM_ARG1, s_divide);
 }
-
-SCM
-scm_divide (SCM x, SCM y)
-{
-  return do_divide (x, y, 0);
-}
-
-static SCM scm_divide2real (SCM x, SCM y)
-{
-  return do_divide (x, y, 1);
-}
 #undef FUNC_NAME
 
 
@@ -9617,12 +9705,11 @@ log_of_fraction (SCM n, SCM d)
                            log_of_exact_integer (d)));
   else if (scm_is_false (scm_negative_p (n)))
     return scm_from_double
-      (log1p (scm_to_double (scm_divide2real (scm_difference (n, d), d))));
+      (log1p (scm_i_divide2double (scm_difference (n, d), d)));
   else
     return scm_c_make_rectangular
-      (log1p (scm_to_double (scm_divide2real
-                            (scm_difference (scm_abs (n), d),
-                             d))),
+      (log1p (scm_i_divide2double (scm_difference (scm_abs (n), d),
+                                   d)),
        M_PI);
 }
 
@@ -9890,6 +9977,16 @@ scm_init_numbers ()
 #endif
 
   exactly_one_half = scm_divide (SCM_INUM1, SCM_I_MAKINUM (2));
+
+  {
+    /* Set scm_i_divide2double_lo2b to (2 b^p - 1) */
+    mpz_init_set_ui (scm_i_divide2double_lo2b, 1);
+    mpz_mul_2exp (scm_i_divide2double_lo2b,
+                  scm_i_divide2double_lo2b,
+                  DBL_MANT_DIG + 1); /* 2 b^p */
+    mpz_sub_ui (scm_i_divide2double_lo2b, scm_i_divide2double_lo2b, 1);
+  }
+
 #include "libguile/numbers.x"
 }
 
-- 
1.7.10.4

>From 8d3f8e1b1fcc1abac1967bb9cf9643b63f3dea6d Mon Sep 17 00:00:00 2001
From: Mark H Weaver <address@hidden>
Date: Tue, 5 Mar 2013 05:45:37 -0500
Subject: [PATCH 09/10] Support reading negative exponents larger than
 SCM_MAXEXP.

* libguile/numbers.c (mem2decimal_from_point): Accept negative exponents
  larger than SCM_MAXEXP that produce subnormals.
---
 libguile/numbers.c |    2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/libguile/numbers.c b/libguile/numbers.c
index 601b2a6..c19b380 100644
--- a/libguile/numbers.c
+++ b/libguile/numbers.c
@@ -5929,7 +5929,7 @@ mem2decimal_from_point (SCM result, SCM mem,
                break;
            }
 
-         if (exponent > SCM_MAXEXP)
+         if (exponent > ((sign == 1) ? SCM_MAXEXP : SCM_MAXEXP + DBL_DIG + 1))
            {
              size_t exp_len = idx - start;
              SCM exp_string = scm_i_substring_copy (mem, start, start + 
exp_len);
-- 
1.7.10.4

>From 9fbdffa44db46dfabb940949f73b142a6a033b0d Mon Sep 17 00:00:00 2001
From: Mark H Weaver <address@hidden>
Date: Tue, 5 Mar 2013 05:47:56 -0500
Subject: [PATCH 10/10] Reimplement idbl2str number printer.

Fixes <http://bugs.gnu.org/13757>.

* libguile/numbers.c (idbl2str): Reimplement.
---
 libguile/numbers.c |  339 +++++++++++++++++++++++++++-------------------------
 1 file changed, 179 insertions(+), 160 deletions(-)

diff --git a/libguile/numbers.c b/libguile/numbers.c
index c19b380..6b3cb16 100644
--- a/libguile/numbers.c
+++ b/libguile/numbers.c
@@ -5267,185 +5267,196 @@ void init_fx_radix(double *fx_list, int radix)
 /* use this array as a way to generate a single digit */
 static const char number_chars[] = "0123456789abcdefghijklmnopqrstuvwxyz";
 
+static mpz_t dbl_minimum_normal_mantissa;
+
 static size_t
-idbl2str (double f, char *a, int radix)
-{
-   int efmt, dpt, d, i, wp;
-   double *fx;
-#ifdef DBL_MIN_10_EXP
-   double f_cpy;
-   int exp_cpy;
-#endif /* DBL_MIN_10_EXP */
-   size_t ch = 0;
-   int exp = 0;
-
-   if(radix < 2 || 
-      radix > SCM_MAX_DBL_RADIX)
-   {
-      /* revert to existing behavior */
-      radix = 10;
-   }
+idbl2str (double dbl, char *a, int radix)
+{
+  int ch = 0;
 
-   wp = scm_dblprec[radix-2];
-   fx = fx_per_radix[radix-2];
+  if (radix < 2 || radix > SCM_MAX_DBL_RADIX)
+    /* revert to existing behavior */
+    radix = 10;
 
-  if (f == 0.0)
+  if (isinf (dbl))
     {
-#ifdef HAVE_COPYSIGN
-      double sgn = copysign (1.0, f);
-
-      if (sgn < 0.0)
-       a[ch++] = '-';
-#endif
-      goto zero;       /*{a[0]='0'; a[1]='.'; a[2]='0'; return 3;} */
+      strcpy (a, (dbl > 0.0) ? "+inf.0" : "-inf.0");
+      return 6;
     }
-
-  if (isinf (f))
+  else if (dbl > 0.0)
+    ;
+  else if (dbl < 0.0)
     {
-      if (f < 0)
-       strcpy (a, "-inf.0");
-      else
-       strcpy (a, "+inf.0");
-      return ch+6;
+      dbl = -dbl;
+      a[ch++] = '-';
     }
-  else if (isnan (f))
+  else if (dbl == 0.0)
     {
-      strcpy (a, "+nan.0");
-      return ch+6;
+      if (!double_is_non_negative_zero (dbl))
+        a[ch++] = '-';
+      strcpy (a + ch, "0.0");
+      return ch + 3;
     }
-
-  if (f < 0.0)
+  else if (isnan (dbl))
     {
-      f = -f;
-      a[ch++] = '-';
+      strcpy (a, "+nan.0");
+      return 6;
     }
 
-#ifdef DBL_MIN_10_EXP  /* Prevent unnormalized values, as from 
-                         make-uniform-vector, from causing infinite loops. */
-  /* just do the checking...if it passes, we do the conversion for our
-     radix again below */
-  f_cpy = f;
-  exp_cpy = exp;
+  /* Algorithm taken from "Printing Floating-Point Numbers Quickly and
+     Accurately" by Robert G. Burger and R. Kent Dybvig */
+  {
+    int e, k;
+    mpz_t f, r, s, mplus, mminus, hi, digit;
+    int f_is_even, f_is_odd;
+    int show_exp = 0;
+
+    mpz_inits (f, r, s, mplus, mminus, hi, digit, NULL);
+    mpz_set_d (f, ldexp (frexp (dbl, &e), DBL_MANT_DIG));
+    if (e < DBL_MIN_EXP)
+      {
+        mpz_tdiv_q_2exp (f, f, DBL_MIN_EXP - e);
+        e = DBL_MIN_EXP;
+      }
+    e -= DBL_MANT_DIG;
 
-  while (f_cpy < 1.0)
-    {
-      f_cpy *= 10.0;
-      if (exp_cpy-- < DBL_MIN_10_EXP)
-       {
-         a[ch++] = '#';
-         a[ch++] = '.';
-         a[ch++] = '#';
-         return ch;
-       }
-    }
-  while (f_cpy > 10.0)
-    {
-      f_cpy *= 0.10;
-      if (exp_cpy++ > DBL_MAX_10_EXP)
-       {
-         a[ch++] = '#';
-         a[ch++] = '.';
-         a[ch++] = '#';
-         return ch;
-       }
-    }
-#endif
+    f_is_even = !mpz_odd_p (f);
+    f_is_odd = !f_is_even;
 
-  while (f < 1.0)
-    {
-      f *= radix;
-      exp--;
-    }
-  while (f > radix)
-    {
-      f /= radix;
-      exp++;
-    }
+    /* Initialize r, s, mplus, and mminus according
+       to Table 1 from the paper. */
+    if (e < 0)
+      {
+        mpz_set_ui (mminus, 1);
+        if (mpz_cmp (f, dbl_minimum_normal_mantissa) != 0
+            || e == DBL_MIN_EXP - DBL_MANT_DIG)
+          {
+            mpz_set_ui (mplus, 1);
+            mpz_mul_2exp (r, f, 1);
+            mpz_mul_2exp (s, mminus, 1 - e);
+          }
+        else
+          {
+            mpz_set_ui (mplus, 2);
+            mpz_mul_2exp (r, f, 2);
+            mpz_mul_2exp (s, mminus, 2 - e);
+          }
+      }
+    else
+      {
+        mpz_set_ui (mminus, 1);
+        mpz_mul_2exp (mminus, mminus, e);
+        if (mpz_cmp (f, dbl_minimum_normal_mantissa) != 0)
+          {
+            mpz_set (mplus, mminus);
+            mpz_mul_2exp (r, f, 1 + e);
+            mpz_set_ui (s, 2);
+          }
+        else
+          {
+            mpz_mul_2exp (mplus, mminus, 1);
+            mpz_mul_2exp (r, f, 2 + e);
+            mpz_set_ui (s, 4);
+          }
+      }
 
-  if (f + fx[wp] >= radix)
+    /* Find the smallest k such that:
+         (r + mplus) / s <  radix^k  (if f is even)
+         (r + mplus) / s <= radix^k  (if f is odd) */
     {
-      f = 1.0;
-      exp++;
-    }
- zero:
-#ifdef ENGNOT 
-  /* adding 9999 makes this equivalent to abs(x) % 3 */
-  dpt = (exp + 9999) % 3;
-  exp -= dpt++;
-  efmt = 1;
-#else
-  efmt = (exp < -3) || (exp > wp + 2);
-  if (!efmt)
-    {
-      if (exp < 0)
-       {
-         a[ch++] = '0';
-         a[ch++] = '.';
-         dpt = exp;
-         while (++dpt)
-           a[ch++] = '0';
-       }
-      else
-       dpt = exp + 1;
+      /* IMPROVE-ME: Make an initial guess to speed this up */
+      mpz_add (hi, r, mplus);
+      k = 0;
+      while (mpz_cmp (hi, s) >= f_is_odd)
+        {
+          mpz_mul_ui (s, s, radix);
+          k++;
+        }
+      if (k == 0)
+        {
+          mpz_mul_ui (hi, hi, radix);
+          while (mpz_cmp (hi, s) < f_is_odd)
+            {
+              mpz_mul_ui (r, r, radix);
+              mpz_mul_ui (mplus, mplus, radix);
+              mpz_mul_ui (mminus, mminus, radix);
+              mpz_mul_ui (hi, hi, radix);
+              k--;
+            }
+        }
     }
-  else
-    dpt = 1;
-#endif
 
-  do
-    {
-      d = f;
-      f -= d;
-      a[ch++] = number_chars[d];
-      if (f < fx[wp])
-       break;
-      if (f + fx[wp] >= 1.0)
-       {
-          a[ch - 1] = number_chars[d+1]; 
-         break;
-       }
-      f *= radix;
-      if (!(--dpt))
-       a[ch++] = '.';
-    }
-  while (wp--);
+    if (k >= 8 || k <= -3)
+      {
+        /* Use scientific notation */
+        show_exp = k - 1;
+        k = 1;
+      }
+    else if (k <= 0)
+      {
+        int i;
 
-  if (dpt > 0)
-    {
-#ifndef ENGNOT
-      if ((dpt > 4) && (exp > 6))
-       {
-         d = (a[0] == '-' ? 2 : 1);
-         for (i = ch++; i > d; i--)
-           a[i] = a[i - 1];
-         a[d] = '.';
-         efmt = 1;
-       }
-      else
-#endif
-       {
-         while (--dpt)
-           a[ch++] = '0';
-         a[ch++] = '.';
-       }
-    }
-  if (a[ch - 1] == '.')
-    a[ch++] = '0';             /* trailing zero */
-  if (efmt && exp)
-    {
-      a[ch++] = 'e';
-      if (exp < 0)
-       {
-         exp = -exp;
-         a[ch++] = '-';
-       }
-      for (i = radix; i <= exp; i *= radix);
-      for (i /= radix; i; i /= radix)
-       {
-          a[ch++] = number_chars[exp / i];
-         exp %= i;
-       }
-    }
+        /* Print leading zeroes */
+        a[ch++] = '0';
+        a[ch++] = '.';
+        for (i = 0; i > k; i--)
+          a[ch++] = '0';
+      }
+
+    for (;;)
+      {
+        int end_1_p, end_2_p;
+        int d;
+
+        mpz_mul_ui (mplus, mplus, radix);
+        mpz_mul_ui (mminus, mminus, radix);
+        mpz_mul_ui (r, r, radix);
+        mpz_fdiv_qr (digit, r, r, s);
+        d = mpz_get_ui (digit);
+
+        mpz_add (hi, r, mplus);
+        end_1_p = (mpz_cmp (r, mminus) < f_is_even);
+        end_2_p = (mpz_cmp (s, hi) < f_is_even);
+        if (end_1_p || end_2_p)
+          {
+            mpz_mul_2exp (r, r, 1);
+            if (!end_2_p)
+              ;
+            else if (!end_1_p)
+              d++;
+            else if (mpz_cmp (r, s) >= !(d & 1))
+              d++;
+            a[ch++] = number_chars[d];
+            if (--k == 0)
+              a[ch++] = '.';
+            break;
+          }
+        else
+          {
+            a[ch++] = number_chars[d];
+            if (--k == 0)
+              a[ch++] = '.';
+          }
+      }
+
+    if (k > 0)
+      {
+        for (; k > 0; k--)
+          a[ch++] = '0';
+        a[ch++] = '.';
+      }
+
+    if (k == 0)
+      a[ch++] = '0';
+
+    if (show_exp)
+      {
+        a[ch++] = 'e';
+        ch += scm_iint2str (show_exp, radix, a + ch);
+      }
+
+    mpz_clears (f, r, s, mplus, mminus, hi, digit, NULL);
+  }
   return ch;
 }
 
@@ -9987,6 +9998,14 @@ scm_init_numbers ()
     mpz_sub_ui (scm_i_divide2double_lo2b, scm_i_divide2double_lo2b, 1);
   }
 
+  {
+    /* Set dbl_minimum_normal_mantissa to b^{p-1} */
+    mpz_init_set_ui (dbl_minimum_normal_mantissa, 1);
+    mpz_mul_2exp (dbl_minimum_normal_mantissa,
+                  dbl_minimum_normal_mantissa,
+                  DBL_MANT_DIG - 1);
+  }
+
 #include "libguile/numbers.x"
 }
 
-- 
1.7.10.4


reply via email to

[Prev in Thread] Current Thread [Next in Thread]