New algorithm, new wrapper

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

1 files changed, 87 insertions, 1 deletions

diff --git a/log2.c b/log2.c
index f462718f4..3ad9988fe 100644
--- a/log2.c
+++ b/log2.c
@@ -29,6 +29,88 @@ MA 02111-1307, USA. */
 mpfr_t __mpfr_const_log2; /* stored value of log(2) with rnd_mode=GMP_RNDZ */
 int __mpfr_const_log2_prec=0; /* precision of stored value */
 
+
+#define A
+#define A1 1
+#define A2 1
+#undef B
+#define C
+#define C1 2
+#define C2 1
+#define NO_FACTORIAL
+#undef R_IS_RATIONAL
+#define GENERIC mpfr_aux_log2
+#include "generic.c" 
+#undef A
+#undef A1
+#undef A2
+#undef NO_FACTORIAL
+#undef GENERIC
+#undef C
+#undef C1
+#undef C2
+
+int
+#if __STDC__
+mpfr_const_aux_log2(mpfr_ptr mylog, mp_rnd_t rnd_mode)
+#else
+mpfr_const_aux_log2(mylog, rnd_mode) mpfr_ptr 
+mylog; mp_rnd_t rnd_mode;
+#endif
+{
+  int prec;
+  mpfr_t tmp1, tmp2, result,tmp3; 
+  mpz_t cst;
+  int good = 0;
+  int logn;
+  int prec_i_want = PREC(mylog);
+  int prec_x;
+  mpz_init(cst);
+  logn =  (int) ceil(log
+                      ((double) PREC(mylog))
+                      /log(2.0)); 
+  prec_x = prec_i_want + logn;
+  while (!good){
+    prec = (int) ceil(log
+		      ((double) (prec_x))
+		      /log(2.0));  
+    mpfr_init2(tmp1, prec_x);
+    mpfr_init2(result, prec_x);
+    mpfr_init2(tmp2, prec_x);
+    mpfr_init2(tmp3, prec_x);
+    mpz_set_ui(cst, 1);
+    mpfr_aux_log2(tmp1, cst, 4, prec-2);
+    mpfr_div_2exp(tmp1, tmp1, 4,GMP_RNDD);
+    mpfr_mul_ui(tmp1, tmp1, 15,GMP_RNDD);
+
+    mpz_set_ui(cst, 3);
+    mpfr_aux_log2(tmp2, cst, 7, prec-2);
+    mpfr_div_2exp(tmp2, tmp2, 7,GMP_RNDD);
+    mpfr_mul_ui(tmp2, tmp2, 5*3,GMP_RNDD);
+    mpfr_sub(result, tmp1, tmp2, GMP_RNDD);
+
+    mpz_set_ui(cst, 13);
+    mpfr_aux_log2(tmp3, cst, 8, prec-2);
+    mpfr_div_2exp(tmp3, tmp3, 8,GMP_RNDD);
+    mpfr_mul_ui(tmp3, tmp3, 3*13,GMP_RNDD);
+    mpfr_sub(result, result, tmp3, GMP_RNDD);
+
+    mpfr_clear(tmp1);
+    mpfr_clear(tmp2);
+    if (mpfr_can_round(result, prec_x, GMP_RNDD, rnd_mode, prec_i_want)){
+      mpfr_set(mylog, result, rnd_mode);
+      mpfr_clear(result);
+      good = 1;
+    } else
+      {
+	mpfr_clear(result);
+	prec_x += logn;
+      }
+  }
+  mpz_clear(cst);
+  return 0;
+}
+
 /* set x to log(2) rounded to precision PREC(x) with direction rnd_mode 
 
    use formula log(2) = sum(1/k/2^k, k=1..infinity)
@@ -63,6 +145,7 @@ mpfr_const_log2(x, rnd_mode) mpfr_ptr x; mp_rnd_t rnd_mode;
       }
   }
 
+  if (precx < 30000){
   /* need to recompute */
   N=2;
   do {
@@ -92,7 +175,10 @@ mpfr_const_log2(x, rnd_mode) mpfr_ptr x; mp_rnd_t rnd_mode;
 #endif
   mpfr_set_z(x, s, rnd_mode);
   EXP(x) -= N;
-
+  } else
+    {
+      mpfr_const_aux_log2(x, rnd_mode);
+    }
   /* stored computed value */
   if (__mpfr_const_log2_prec==0) mpfr_init2(__mpfr_const_log2, precx);
   else mpfr_set_prec(__mpfr_const_log2, precx);