improved the computation of g = ceil((e-1)/log_2(beta)), using two tables

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

1 files changed, 140 insertions, 38 deletions

diff --git a/get_str.c b/get_str.c
index 600f48af6..f71d931b3 100644
--- a/get_str.c
+++ b/get_str.c
@@ -34,46 +34,135 @@ static double _mpfr_ceil _PROTO ((double));
 static long mpn_exp _PROTO ((mp_limb_t *, mp_exp_t *, int, mp_exp_t, size_t));
 static int mpfr_get_str_aux _PROTO ((char *, mp_exp_t *, mp_limb_t *,
 		       mp_size_t, mp_exp_t, long, int, size_t, mp_rnd_t));
+static mp_exp_t mpfr_get_str_compute_g _PROTO ((int, mp_exp_t));
 
 static char num_to_text[] = "0123456789abcdefghijklmnopqrstuvwxyz";
 
-/* log_b2[b-2] = floor(log(2)/log(b)) for 2 <= b <= 36 */
-static const double log_b2[35] = 
-  {0.10000000000000000000E1,
-   0.63092975357145743709E0,
-   0.50000000000000000000E0,
-   0.43067655807339305067E0,
-   0.38685280723454158687E0,
-   0.35620718710802217651E0,
-   0.33333333333333333333E0,
-   0.31546487678572871854E0,
-   0.30102999566398119521E0,
-   0.28906482631788785926E0,
-   0.27894294565112984319E0,
-   0.27023815442731974129E0,
-   0.26264953503719354797E0,
-   0.25595802480981548938E0,
-   0.25000000000000000000E0,
-   0.24465054211822603038E0,
-   0.23981246656813144473E0,
-   0.23540891336663823644E0,
-   0.23137821315975917426E0,
-   0.22767024869695299798E0,
-   0.22424382421757543947E0,
-   0.22106472945750374614E0,
-   0.21810429198553155922E0,
-   0.21533827903669652533E0,
-   0.21274605355336315360E0,
-   0.21030991785715247903E0,
-   0.20801459767650945759E0,
-   0.20584683246043445730E0,
-   0.20379504709050619003E0,
-   0.20184908658209985071E0,
-   0.19999999999999999999E0,
-   0.19823986317056053160E0,
-   0.19656163223282260771E0,
-   0.19495902189378630772E0,
-   0.19342640361727079343E0};
+/* 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>
+#include "gmp.h"
+#include "mpfr.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_round_prec (l, GMP_RNDU, 53);
+      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} */
@@ -402,6 +491,19 @@ mpfr_get_str_aux (char *str, mp_exp_t *exp, mp_limb_t *r, mp_size_t n,
   return (dir);
 }
 
+/* returns ceil(e/log_2(beta)) */
+static mp_exp_t
+mpfr_get_str_compute_g (int beta, mp_exp_t e)
+{
+  double g0, g1;
+  mp_exp_t g;
+  
+  g0 = (double) e * log_b2[beta - 2];
+  g1 = (double) e * log_b2_low[beta - 2];
+  g = (mp_exp_t) _mpfr_ceil (g0);
+  g0 -= (double) g;
+  return g + _mpfr_ceil (g0 + g1);
+}
 
 /* prints the mantissa of x in the string s, and writes the corresponding
    exponent in e.
@@ -571,7 +673,7 @@ mpfr_get_str (char *s, mp_exp_t *e, int b, size_t m, mpfr_srcptr x, mp_rnd_t rnd
       return (s0);
    }
 
-  g = (size_t) _mpfr_ceil (((double) MPFR_EXP(x) - 1.0) * log_b2[b - 2]);
+  g = mpfr_get_str_compute_g (b, MPFR_EXP(x) - 1);
   exact = 1;
   prec = (mp_exp_t) _mpfr_ceil ((double) m / log_b2[b-2]) + 1;
   exp = ((mp_exp_t) m < g) ? g - (mp_exp_t) m : (mp_exp_t) m - g;