[Guile-commits] GNU Guile branch, stable-2.0, updated. v2.0.7-204-g1ea3762

[Mark H Weaver]

This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "GNU Guile".

http://git.savannah.gnu.org/cgit/guile.git/commit/?id=1ea37620c2c1794f7685b312d2530676a078ada7

The branch, stable-2.0 has been updated
       via  1ea37620c2c1794f7685b312d2530676a078ada7 (commit)
      from  982377849029f2840ebb105cda49390fecca4fe4 (commit)

Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.

- Log -----------------------------------------------------------------
commit 1ea37620c2c1794f7685b312d2530676a078ada7
Author: Mark H Weaver <address@hidden>
Date:   Tue Mar 5 05:47:56 2013 -0500

    Reimplement idbl2str number printer.
    
    Fixes <http://bugs.gnu.org/13757>.
    
    * libguile/numbers.c (idbl2str): Reimplement.
      (mem2decimal_from_point): Accept negative exponents larger than
      SCM_MAXEXP that produce subnormals.
      (SCM_MAX_DBL_PREC): Removed preprocessor macro.
      (scm_dblprec, fx_per_radix): Removed static variables.
      (init_dblprec, init_fx_radix): Removed static functions.
      (scm_init_numbers): Remove initialization code for 'scm_dblprec'
      and 'fx_per_radix'.
    
    * test-suite/tests/numbers.test ("number->string"): Restore tests that
      previously failed.  Remove comments about problems in the number
      printer that are now fixed.

-----------------------------------------------------------------------

Summary of changes:
 libguile/numbers.c            |  391 +++++++++++++++++++----------------------
 test-suite/tests/numbers.test |   97 ++++++++---
 2 files changed, 254 insertions(+), 234 deletions(-)

diff --git a/libguile/numbers.c b/libguile/numbers.c
index f632733..c641e3f 100644
--- a/libguile/numbers.c
+++ b/libguile/numbers.c
@@ -5250,229 +5250,201 @@ SCM_DEFINE (scm_integer_length, "integer-length", 1, 
0, 0,
 #undef FUNC_NAME
 
 /*** NUMBERS -> STRINGS ***/
-#define SCM_MAX_DBL_PREC  60
 #define SCM_MAX_DBL_RADIX 36
 
-/* this is an array starting with radix 2, and ending with radix 
SCM_MAX_DBL_RADIX */
-static int scm_dblprec[SCM_MAX_DBL_RADIX - 1];
-static double fx_per_radix[SCM_MAX_DBL_RADIX - 1][SCM_MAX_DBL_PREC];
-
-static
-void init_dblprec(int *prec, int radix) {
-   /* determine floating point precision by adding successively
-      smaller increments to 1.0 until it is considered == 1.0 */
-   double f = ((double)1.0)/radix;
-   double fsum = 1.0 + f;
-
-   *prec = 0;
-   while (fsum != 1.0)
-   {
-      if (++(*prec) > SCM_MAX_DBL_PREC)
-         fsum = 1.0;
-      else
-      {
-         f /= radix;
-         fsum = f + 1.0;
-      }
-   }
-   (*prec) -= 1;
-}
-
-static
-void init_fx_radix(double *fx_list, int radix)
-{
-  /* initialize a per-radix list of tolerances.  When added
-     to a number < 1.0, we can determine if we should raund
-     up and quit converting a number to a string. */
-   int i;
-   fx_list[0] = 0.0;
-   fx_list[1] = 0.5;
-   for( i = 2 ; i < SCM_MAX_DBL_PREC; ++i ) 
-     fx_list[i] = (fx_list[i-1] / 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)
+idbl2str (double dbl, 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;
-   }
+  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;
 }
 
@@ -5956,7 +5928,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);
@@ -9993,8 +9965,6 @@ SCM_PRIMITIVE_GENERIC (scm_sqrt, "sqrt", 1, 0, 0,
 void
 scm_init_numbers ()
 {
-  int i;
-
   if (scm_install_gmp_memory_functions)
     mp_set_memory_functions (custom_gmp_malloc,
                              custom_gmp_realloc,
@@ -10016,17 +9986,6 @@ scm_init_numbers ()
   flo0 = scm_from_double (0.0);
   flo_log10e = scm_from_double (M_LOG10E);
 
-  /* determine floating point precision */
-  for (i=2; i <= SCM_MAX_DBL_RADIX; ++i)
-    {
-      init_dblprec(&scm_dblprec[i-2],i);
-      init_fx_radix(fx_per_radix[i-2],i);
-    }
-#ifdef DBL_DIG
-  /* hard code precision for base 10 if the preprocessor tells us to... */
-  scm_dblprec[10-2] = (DBL_DIG > 20) ? 20 : DBL_DIG;
-#endif
-
   exactly_one_half = scm_divide (SCM_INUM1, SCM_I_MAKINUM (2));
 
   {
@@ -10038,6 +9997,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"
 }
 
diff --git a/test-suite/tests/numbers.test b/test-suite/tests/numbers.test
index 5a77e93..8f01633 100644
--- a/test-suite/tests/numbers.test
+++ b/test-suite/tests/numbers.test
@@ -1383,21 +1383,39 @@
          (lambda (n radix)
            (string->number (number->string n radix) radix))))
 
+    (define (test num)
+      (pass-if-equal (list num 'pos)
+          num
+        (num->str->num num 10))
+      (pass-if-equal (list num 'neg)
+          (- num)
+        (num->str->num (- num) 10)))
+
     (pass-if (documented? number->string))
     (pass-if (string=? (number->string 0) "0"))
     (pass-if (string=? (number->string 171) "171"))
     (pass-if (= (+ fixnum-max 1) (num->str->num (+ fixnum-max 1) 10)))
     (pass-if (= (- fixnum-min 1) (num->str->num (- fixnum-min 1) 10)))
-    (pass-if (= (inf) (num->str->num (inf) 10)))
-    (pass-if (= 1.3 (num->str->num 1.3 10)))
 
-    ;; XXX - some results depend on whether Guile is compiled optimzed
-    ;; or not.  It is clearly undesirable to have number->string to be
-    ;; influenced by this.
+    (test (inf))
+    (test 1.3)
+    (test (acos -1))  ; pi
+    (test (exp 1))    ; e
+    (test (/ 3.0))
+    (test (/ 7.0))
+    (test 2.2250738585072011e-308)
+    (test 2.2250738585072012e-308)
+
+    ;; Largest finite inexact
+    (test (* (- (expt 2.0 dbl-mant-dig) 1)
+             (expt 2.0 (- dbl-max-exp dbl-mant-dig))))
+
+    (pass-if (string=? "0.0" (number->string 0.0)))
+    (pass-if (or (eqv? 0.0 -0.0)
+                 (string=? "-0.0" (number->string -0.0))))
 
     (pass-if (string=? (number->string 35.25 36) "z.9"))
-    (pass-if (or (string=? (number->string 0.25 2) "0.01")
-                (string=? (number->string 0.25 2) "0.010")))
+    (pass-if (string=? (number->string 0.25 2) "0.01"))
     (pass-if (string=? (number->string 255.0625 16) "ff.1"))
     (pass-if (string=? (number->string (/ 1 3) 3) "1/10"))
 
@@ -1411,26 +1429,61 @@
     (pass-if (string=? (number->string 35 36) "z"))
     (pass-if (= (num->str->num 35 36) 35))
 
+    (with-test-prefix "powers of radix"
+      (for-each
+       (lambda (radix)
+         (for-each (lambda (k)
+                     (let ((val (exact->inexact (expt radix k)))
+                           (str (if (<= -3 k 6)
+                                    (assoc-ref '((-3 . "0.001")
+                                                 (-2 . "0.01")
+                                                 (-1 . "0.1")
+                                                 ( 0 . "1.0")
+                                                 ( 1 . "10.0")
+                                                 ( 2 . "100.0")
+                                                 ( 3 . "1000.0")
+                                                 ( 4 . "10000.0")
+                                                 ( 5 . "100000.0")
+                                                 ( 6 . "1000000.0"))
+                                               k)
+                                    (string-append "1.0e"
+                                                   (number->string k radix)))))
+                       (pass-if-equal (list radix k 'pos)
+                           str
+                         (number->string val radix))
+                       (pass-if-equal (list radix k 'neg)
+                           (string-append "-" str)
+                         (number->string (- val) radix))))
+                   (iota 41 -20)))
+       (iota 35 2)))
+
+    (with-test-prefix "multiples of smallest inexact"
+      (for-each (lambda (k)
+                  (let ((val (* k (expt 2.0 (- dbl-min-exp dbl-mant-dig)))))
+                    (test val)))
+                (iota 40 1)))
+
+    (with-test-prefix "one plus multiples of epsilon"
+      (for-each (lambda (k)
+                  (let ((val (+ 1.0 (* k dbl-epsilon))))
+                    (test val)))
+                (iota 40 1)))
+
+    (with-test-prefix "one minus multiples of 1/2 epsilon"
+      (for-each (lambda (k)
+                  (let ((val (- 1.0 (* k 1/2 dbl-epsilon))))
+                    (test val)))
+                (iota 40 1)))
+
     ;; Before Guile 2.0.1, even in the presence of a #e forced exactness
     ;; specifier, negative exponents were applied inexactly and then
     ;; later coerced to exact, yielding an incorrect fraction.
     (pass-if (eqv? (string->number "#e1e-10") 1/10000000000))
 
-    ;; Numeric conversion from decimal is not precise, in its current
-    ;; implementation, so 11.333... and 1.324... can't be expected to
-    ;; reliably come out to precise values.  These tests did actually work
-    ;; for a while, but something in gcc changed, affecting the conversion
-    ;; code.
-    ;;
-    ;; (pass-if (or (string=? (number->string 11.33333333333333333 12)
-    ;;                        "B.4")
-    ;;              (string=? (number->string 11.33333333333333333 12)
-    ;;                        "B.400000000000009")))
-    ;; (pass-if (or (string=? (number->string 1.324e44 16)
-    ;;                        "5.EFE0A14FAFEe24")
-    ;;              (string=? (number->string 1.324e44 16)
-    ;;                        "5.EFE0A14FAFDF8e24")))
-    ))
+    (pass-if (string=? (number->string 11.33333333333333333 12)
+                       "b.4"))
+    (pass-if (string=? (number->string 1.324e44 16)
+                       "5.efe0a14fafdf8e24"))))
   
 ;;;
 ;;; string->number


hooks/post-receive
-- 
GNU Guile