Replaced some computations using the type double by computations using

MPFR, for mpfr_get_str, allowing it to work with an x86 processor set
up in single-precision mode.


git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@4333 280ebfd0-de03-0410-8827-d642c229c3f4

1 files changed, 38 insertions, 157 deletions

diff --git a/get_str.c b/get_str.c
index 5735e9abf..2e4fed72b 100644
--- a/get_str.c
+++ b/get_str.c
@@ -21,7 +21,7 @@ the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
 MA 02110-1301, USA. */
 
 #include <string.h> /* For strlen */
-#include <limits.h> /* For LONG_MAX and LONG_MIN */
+#include <limits.h> /* For CHAR_BIT, LONG_MAX and LONG_MIN */
 
 #define MPFR_NEED_LONGLONG_H
 #include "mpfr-impl.h"
@@ -29,134 +29,9 @@ MA 02110-1301, USA. */
 static double mpfr_ceil_double (double);
 static int mpfr_get_str_aux (char *const, mp_exp_t *const, mp_limb_t *const,
                        mp_size_t, mp_exp_t, long, int, size_t, mp_rnd_t);
-static mp_exp_t mpfr_get_str_compute_g (int, mp_exp_t);
 
 static const char num_to_text[] = "0123456789abcdefghijklmnopqrstuvwxyz";
 
-/* for 2 <= b <= 36, log_b2[b-2] + log_b2_low[b-2] is a 76-bit upper
-   approximation of log(2)/log(b), with log_b2[b-2] having 23 significative
-   bits only. These approximations were computed with the following program.
-
-#include <stdio.h>
-
-double log_b2[35], log_b2_low[35];
-
-main()
-{
-  int beta;
-  mpfr_t l, l0;
-
-  for (beta=2;beta<=36;beta++)
-    {
-      mpfr_init2 (l, 77);
-      mpfr_set_ui (l, beta, GMP_RNDD);
-      mpfr_log2 (l, l, GMP_RNDD);
-      mpfr_ui_div (l, 1, l, GMP_RNDU);
-      mpfr_init2 (l0, 23);
-      mpfr_set (l0, l, GMP_RNDD);
-      mpfr_sub (l, l, l0, GMP_RNDU);
-      mpfr_prec_round (l, 53, GMP_RNDU);
-      log_b2[beta-2] = mpfr_get_d (l0, GMP_RNDU);
-      log_b2_low[beta-2] = mpfr_get_d (l, GMP_RNDU);
-      mpfr_clear (l0);
-      mpfr_clear (l);
-    }
-
-  printf ("static const double log_b2[35] = {");
-  for (beta=2;beta<=36;beta++)
-    {
-      printf ("\n%1.20e", log_b2[beta-2]);
-      if (beta < 36) printf (",");
-    }
-  printf ("\n};\n");
-
-  printf ("static const double log_b2_low[35] = {");
-  for (beta=2;beta<=36;beta++)
-    {
-      printf ("\n%1.20e", log_b2_low[beta-2]);
-      if (beta < 36) printf (",");
-    }
-  printf ("\n};\n");
-}
-*/
-
-
-static const double log_b2[35] = {
-  1.00000000000000000000e+00,
-  6.30929708480834960938e-01,
-  5.00000000000000000000e-01,
-  4.30676519870758056641e-01,
-  3.86852800846099853516e-01,
-  3.56207132339477539062e-01,
-  3.33333313465118408203e-01,
-  3.15464854240417480469e-01,
-  3.01029980182647705078e-01,
-  2.89064824581146240234e-01,
-  2.78942942619323730469e-01,
-  2.70238101482391357422e-01,
-  2.62649476528167724609e-01,
-  2.55958020687103271484e-01,
-  2.50000000000000000000e-01,
-  2.44650512933731079102e-01,
-  2.39812463521957397461e-01,
-  2.35408902168273925781e-01,
-  2.31378197669982910156e-01,
-  2.27670222520828247070e-01,
-  2.24243819713592529297e-01,
-  2.21064716577529907227e-01,
-  2.18104273080825805664e-01,
-  2.15338259935379028320e-01,
-  2.12746024131774902344e-01,
-  2.10309892892837524414e-01,
-  2.08014577627182006836e-01,
-  2.05846816301345825195e-01,
-  2.03795045614242553711e-01,
-  2.01849073171615600586e-01,
-  1.99999988079071044922e-01,
-  1.98239862918853759766e-01,
-  1.96561604738235473633e-01,
-  1.94959014654159545898e-01,
-  1.93426400423049926758e-01
-};
-
-static const double log_b2_low[35] = {
-  0.00000000000000000000e+00,
-  4.50906224761620348192e-08,
-  0.00000000000000000000e+00,
-  3.82026349940294905572e-08,
-  6.38844173335462308442e-09,
-  5.47685446374516835508e-08,
-  1.98682149251302105391e-08,
-  2.25453112380810174096e-08,
-  1.54813334901356166450e-08,
-  1.73674161903183898500e-09,
-  3.03180611272229558617e-09,
-  5.29449283838722401896e-08,
-  5.85090258233695360801e-08,
-  4.12271221790330610183e-09,
-  0.00000000000000000000e+00,
-  2.91844949512882013763e-08,
-  3.04617404727474892263e-09,
-  1.11983643106657188415e-08,
-  1.54897762641074502508e-08,
-  2.61761247509117662805e-08,
-  4.50398291018069050027e-09,
-  1.28799738389232369510e-08,
-  1.89047057535652783545e-08,
-  1.91013174970147452786e-08,
-  2.94215882512624420131e-08,
-  2.49643149546191168760e-08,
-  2.00493274507626165734e-08,
-  1.61590886321114899524e-08,
-  1.47626363632181082532e-09,
-  1.34104842501325808254e-08,
-  1.19209289550781256617e-08,
-  2.51706771840338761560e-10,
-  2.74945871340834855649e-08,
-  7.23962676182708790191e-09,
-  3.19422086667731154221e-09
-};
-
 /* copy most important limbs of {op, n2} in {rp, n1} */
 /* if n1 > n2 put 0 in low limbs of {rp, n1} */
 #define MPN_COPY2(rp, n1, op, n2) \
@@ -170,17 +45,6 @@ static const double log_b2_low[35] = {
       MPN_ZERO ((rp), (n1) - (n2)); \
     }
 
-static double
-mpfr_ceil_double (double x)
-{
-  double y;
-  /* Note: this function should be rewritten to avoid the possible
-     overflow. */
-  MPFR_ASSERTN(x >= (double) LONG_MIN && x <= (double) LONG_MAX);
-  y = (double) (long int) x;
-  return ((y < x) ? y + 1.0 : y);
-}
-
 #define MPFR_ROUND_FAILED 3
 
 /* Input: an approximation r*2^f of an real Y, with |r*2^f-Y| <= 2^(e+f).
@@ -369,26 +233,44 @@ mpfr_get_str_aux (char *const str, mp_exp_t *const exp, mp_limb_t *const r,
   return dir;
 }
 
-/* returns ceil(e/log_2(beta)) */
+/* mpfr_l2b[b-2][0] is a 23-bit upper approximation to log(b)/log(2),
+   mpfr_l2b[b-2][1] is a 76-bit upper approximation to log(2)/log(b),
+   if not null. The array is initialized with null pointers. */
+mpfr_ptr MPFR_THREAD_ATTR mpfr_l2b[BASE_MAX-1][2] = { NULL };
+
+/* returns ceil(e * log2(b)^((-1)^i)), or ... + 1 */
 static mp_exp_t
-mpfr_get_str_compute_g (int beta, mp_exp_t e)
+ceil_mul (mp_exp_t e, int beta, int i)
 {
-  double g0, g1, de;
-  mp_exp_t g;
-
-  de = (double) e;
-  g0 = de * log_b2[beta - 2];
-  g1 = de * log_b2_low[beta - 2];
-  if (de > 9007199254740992.0 || de < -9007199254740992.0)
-    /* can happen on 64-bit machines */
+  mpfr_ptr p;
+  mpfr_t t;
+  mp_exp_t r;
+
+  p = mpfr_l2b[beta-2][i];
+  if (p == NULL)
     {
-      mp_exp_t low_e = e - (mp_exp_t) de;
-      g1 += (double) low_e * log_b2[beta - 2];
+      mpfr_l2b[beta-2][i] = p =
+        (mpfr_ptr) (*__gmp_allocate_func) (sizeof (mpfr_t));
+      if (i == 0)
+        {
+          mpfr_init2 (p, 23);
+          mpfr_set_ui (p, beta, GMP_RNDU);
+          mpfr_log2 (p, p, GMP_RNDU);
+        }
+      else
+        {
+          mpfr_init2 (p, 77);
+          mpfr_set_ui (p, beta, GMP_RNDD);
+          mpfr_log2 (p, p, GMP_RNDD);
+          mpfr_ui_div (p, 1, p, GMP_RNDU);
+        }
     }
-  g = (mp_exp_t) mpfr_ceil_double (g0);
-  g0 -= (double) g;
-  g += (mp_exp_t) mpfr_ceil_double (g0 + g1);
-  return g;
+  mpfr_init2 (t, sizeof (mp_exp_t) * CHAR_BIT);
+  mpfr_set_exp_t (t, e, GMP_RNDU);
+  mpfr_mul (t, t, p, GMP_RNDU);
+  r = mpfr_get_exp_t (t, GMP_RNDU);
+  mpfr_clear (t);
+  return r;
 }
 
 /* prints the mantissa of x in the string s, and writes the corresponding
@@ -442,8 +324,7 @@ mpfr_get_str (char *s, mp_exp_t *e, int b, size_t m, mpfr_srcptr x, mp_rnd_t rnd
          the first base-b digit contains only one bit, so we get
          1 + ceil((n-1)/k) = 2 + floor((n-2)/k) instead.
       */
-      m = 1 + mpfr_get_str_compute_g (b, (IS_POW2(b)) ? MPFR_PREC(x) - 1
-                                      : MPFR_PREC(x));
+      m = 1 + ceil_mul (IS_POW2(b) ? MPFR_PREC(x) - 1: MPFR_PREC(x), b, 1);
       if (m < 2)
         m = 2;
     }
@@ -564,9 +445,9 @@ mpfr_get_str (char *s, mp_exp_t *e, int b, size_t m, mpfr_srcptr x, mp_rnd_t rnd
   if (neg)
     rnd = MPFR_INVERT_RND(rnd);
 
-  g = mpfr_get_str_compute_g (b, MPFR_GET_EXP (x) - 1);
+  g = ceil_mul (MPFR_GET_EXP (x) - 1, b, 1);
   exact = 1;
-  prec = (mp_exp_t) mpfr_ceil_double ((double) m / log_b2[b-2]) + 1;
+  prec = ceil_mul (m, b, 0) + 1;
   exp = ((mp_exp_t) m < g) ? g - (mp_exp_t) m : (mp_exp_t) m - g;
   prec += MPFR_INT_CEIL_LOG2 (prec); /* number of guard bits */
   if (exp != 0) /* add maximal exponentiation error */