Macros MPFR_EXP_INVALID (invalid exponent value) and MPFR_EXP_CHECK

added. Code update to use MPFR_GET_EXP and MPFR_SET_EXP instead of
MPFR_EXP to allow more bug detection related to special values.
Macros MPFR_SET_NAN, MPFR_SET_INF, MPFR_SET_ZERO and MPFR_INIT set
the exponent of the number to MPFR_EXP_INVALID if MPFR_EXP_CHECK
is defined. Compile with -DMPFR_EXP_CHECK and make check to see
the potential problems; currently, 40 of 76 tests fail.


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

1 files changed, 127 insertions, 103 deletions

diff --git a/const_pi.c b/const_pi.c
index f50581790..d57c201e6 100644
--- a/const_pi.c
+++ b/const_pi.c
@@ -1,6 +1,6 @@
 /* mpfr_const_pi -- compute Pi
 
-Copyright 1999, 2000, 2001, 2002 Free Software Foundation.
+Copyright 1999, 2000, 2001, 2002, 2003 Free Software Foundation.
 
 This file is part of the MPFR Library.
 
@@ -53,67 +53,70 @@ mpfr_pi_machin3 (mpfr_ptr mylog, mp_rnd_t rnd_mode)
   mpz_t cst;
 
   MPFR_CLEAR_FLAGS(mylog); 
-  logn =  __gmpfr_ceil_log2 ((double) MPFR_PREC(mylog));
+  logn = __gmpfr_ceil_log2 ((double) MPFR_PREC(mylog));
   prec_x = prec_i_want + logn + 5;
   mpz_init(cst);  
-  while (!good){
-  prec = __gmpfr_ceil_log2 ((double) prec_x);
-
-  mpfr_init2(tmp1, prec_x);
-  mpfr_init2(tmp2, prec_x);
-  mpfr_init2(tmp3, prec_x);
-  mpfr_init2(tmp4, prec_x);
-  mpfr_init2(tmp5, prec_x);
-  mpfr_init2(tmp6, prec_x);
-  mpfr_init2(result, prec_x);
-  mpz_set_si(cst, -1);
-
-  mpfr_aux_pi(tmp1, cst, 268*268, prec - 4);
-  mpfr_div_ui(tmp1, tmp1, 268, GMP_RNDD);
-  mpfr_mul_ui(tmp1, tmp1, 61, GMP_RNDD);
-
-  mpfr_aux_pi(tmp2, cst, 343*343, prec - 4);
-  mpfr_div_ui(tmp2, tmp2, 343, GMP_RNDD);
-  mpfr_mul_ui(tmp2, tmp2, 122, GMP_RNDD);
-
-  mpfr_aux_pi(tmp3, cst, 557*557, prec - 4);
-  mpfr_div_ui(tmp3, tmp3, 557, GMP_RNDD);
-  mpfr_mul_ui(tmp3, tmp3, 115, GMP_RNDD);
-
-  mpfr_aux_pi(tmp4, cst, 1068*1068, prec - 4);
-  mpfr_div_ui(tmp4, tmp4, 1068, GMP_RNDD);
-  mpfr_mul_ui(tmp4, tmp4, 32, GMP_RNDD);
-
-  mpfr_aux_pi(tmp5, cst, 3458*3458, prec - 4);
-  mpfr_div_ui(tmp5, tmp5, 3458, GMP_RNDD);
-  mpfr_mul_ui(tmp5, tmp5, 83, GMP_RNDD);
-
-  mpfr_aux_pi(tmp6, cst, 27493*27493, prec - 4);
-  mpfr_div_ui(tmp6, tmp6, 27493, GMP_RNDD);
-  mpfr_mul_ui(tmp6, tmp6, 44, GMP_RNDD);
-
-  mpfr_add(result, tmp1, tmp2, GMP_RNDD);
-  mpfr_add(result, result, tmp3, GMP_RNDD);
-  mpfr_sub(result, result, tmp4, GMP_RNDD);
-  mpfr_add(result, result, tmp5, GMP_RNDD);
-  mpfr_add(result, result, tmp6, GMP_RNDD);
-  mpfr_mul_2ui(result, result, 2, GMP_RNDD);
-  mpfr_clear(tmp1);
-  mpfr_clear(tmp2);
-  mpfr_clear(tmp3);
-  mpfr_clear(tmp4);
-  mpfr_clear(tmp5);
-  mpfr_clear(tmp6);
-  if (mpfr_can_round(result, prec_x - 5, GMP_RNDD, rnd_mode, prec_i_want)){
-    mpfr_set(mylog, result, rnd_mode);
-    mpfr_clear(result);
-    good = 1;
-  } else
+  while (!good)
     {
-      mpfr_clear(result);
-      prec_x += logn;
-      }
-  }
+      prec = __gmpfr_ceil_log2 ((double) prec_x);
+
+      mpfr_init2(tmp1, prec_x);
+      mpfr_init2(tmp2, prec_x);
+      mpfr_init2(tmp3, prec_x);
+      mpfr_init2(tmp4, prec_x);
+      mpfr_init2(tmp5, prec_x);
+      mpfr_init2(tmp6, prec_x);
+      mpfr_init2(result, prec_x);
+      mpz_set_si(cst, -1);
+
+      mpfr_aux_pi(tmp1, cst, 268*268, prec - 4);
+      mpfr_div_ui(tmp1, tmp1, 268, GMP_RNDD);
+      mpfr_mul_ui(tmp1, tmp1, 61, GMP_RNDD);
+
+      mpfr_aux_pi(tmp2, cst, 343*343, prec - 4);
+      mpfr_div_ui(tmp2, tmp2, 343, GMP_RNDD);
+      mpfr_mul_ui(tmp2, tmp2, 122, GMP_RNDD);
+
+      mpfr_aux_pi(tmp3, cst, 557*557, prec - 4);
+      mpfr_div_ui(tmp3, tmp3, 557, GMP_RNDD);
+      mpfr_mul_ui(tmp3, tmp3, 115, GMP_RNDD);
+
+      mpfr_aux_pi(tmp4, cst, 1068*1068, prec - 4);
+      mpfr_div_ui(tmp4, tmp4, 1068, GMP_RNDD);
+      mpfr_mul_ui(tmp4, tmp4, 32, GMP_RNDD);
+
+      mpfr_aux_pi(tmp5, cst, 3458*3458, prec - 4);
+      mpfr_div_ui(tmp5, tmp5, 3458, GMP_RNDD);
+      mpfr_mul_ui(tmp5, tmp5, 83, GMP_RNDD);
+
+      mpfr_aux_pi(tmp6, cst, 27493*27493, prec - 4);
+      mpfr_div_ui(tmp6, tmp6, 27493, GMP_RNDD);
+      mpfr_mul_ui(tmp6, tmp6, 44, GMP_RNDD);
+
+      mpfr_add(result, tmp1, tmp2, GMP_RNDD);
+      mpfr_add(result, result, tmp3, GMP_RNDD);
+      mpfr_sub(result, result, tmp4, GMP_RNDD);
+      mpfr_add(result, result, tmp5, GMP_RNDD);
+      mpfr_add(result, result, tmp6, GMP_RNDD);
+      mpfr_mul_2ui(result, result, 2, GMP_RNDD);
+      mpfr_clear(tmp1);
+      mpfr_clear(tmp2);
+      mpfr_clear(tmp3);
+      mpfr_clear(tmp4);
+      mpfr_clear(tmp5);
+      mpfr_clear(tmp6);
+      if (mpfr_can_round(result, prec_x - 5, 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;    
 }
@@ -168,52 +171,73 @@ mpfr_const_pi (mpfr_ptr x, mp_rnd_t rnd_mode)
       mpfr_set(x, __mpfr_const_pi, rnd_mode); return; 
     }
 
-  if (prec < 20000){
-    /* need to recompute */
-    N=1; 
-    do {
-      oldN = N;
-      N = (prec+3)/4 + __gmpfr_ceil_log2((double) N + 1.0);
-    } while (N != oldN);
-    mpz_init(pi); mpz_init(num); mpz_init(den); mpz_init(d3); mpz_init(d2);
-    mpz_init(tmp);
-    mpz_set_ui(pi, 0);
-    mpz_set_ui(num, 16); /* num(-1) */
-    mpz_set_ui(den, 21); /* den(-1) */
-    mpz_set_si(d3, -2454);
-    mpz_set_ui(d2, 14736);
-  /* invariants: num=120*n^2+151*n+47, den=512*n^4+1024*n^3+712*n^2+194*n+15
-     d3 = 2048*n^3+400*n-6, d2 = 6144*n^2-6144*n+2448
-   */
-    for (n=0; n<N; n++) {
-      /* num(n)-num(n-1) = 240*n+31 */
-      mpz_add_ui(num, num, 240*n+31); /* no overflow up to MPFR_PREC=71M */
-      /* d2(n) - d2(n-1) = 12288*(n-1) */
-      if (n>0) mpz_add_ui(d2, d2, 12288*(n-1));
-      else mpz_sub_ui(d2, d2, 12288);
-      /* d3(n) - d3(n-1) = d2 */
-      mpz_add(d3, d3, d2);
-      /* den(n)-den(n-1) = 2048*n^3 + 400n - 6 = d3 */
-      mpz_add(den, den, d3);
-      mpz_mul_2exp(tmp, num, 4*(N-n));
-      mpz_fdiv_q(tmp, tmp, den);
-      mpz_add(pi, pi, tmp);
-    }
-    mpfr_set_z(x, pi, rnd_mode);
-    mpfr_init2(y, mpfr_get_prec(x));
-    mpz_add_ui(pi, pi, N+1);
-    mpfr_set_z(y, pi, rnd_mode);
-    if (mpfr_cmp(x, y) != 0) {
-      fprintf(stderr, "does not converge\n"); exit(1);
+  if (prec < 20000)
+    {
+      /* need to recompute */
+      N=1; 
+      do
+        {
+          oldN = N;
+          N = (prec+3)/4 + __gmpfr_ceil_log2((double) N + 1.0);
+        }
+      while (N != oldN);
+      mpz_init(pi);
+      mpz_init(num);
+      mpz_init(den);
+      mpz_init(d3);
+      mpz_init(d2);
+      mpz_init(tmp);
+      mpz_set_ui(pi, 0);
+      mpz_set_ui(num, 16); /* num(-1) */
+      mpz_set_ui(den, 21); /* den(-1) */
+      mpz_set_si(d3, -2454);
+      mpz_set_ui(d2, 14736);
+      /* invariants: num=120*n^2+151*n+47,
+         den=512*n^4+1024*n^3+712*n^2+194*n+15
+         d3 = 2048*n^3+400*n-6, d2 = 6144*n^2-6144*n+2448
+      */
+      for (n=0; n<N; n++)
+        {
+          /* num(n)-num(n-1) = 240*n+31 */
+          mpz_add_ui(num, num, 240*n+31); /* no overflow up to MPFR_PREC=71M */
+          /* d2(n) - d2(n-1) = 12288*(n-1) */
+          if (n>0)
+            mpz_add_ui(d2, d2, 12288*(n-1));
+          else
+            mpz_sub_ui(d2, d2, 12288);
+          /* d3(n) - d3(n-1) = d2 */
+          mpz_add(d3, d3, d2);
+          /* den(n)-den(n-1) = 2048*n^3 + 400n - 6 = d3 */
+          mpz_add(den, den, d3);
+          mpz_mul_2exp(tmp, num, 4*(N-n));
+          mpz_fdiv_q(tmp, tmp, den);
+          mpz_add(pi, pi, tmp);
+        }
+      mpfr_set_z(x, pi, rnd_mode);
+      mpfr_init2(y, mpfr_get_prec(x));
+      mpz_add_ui(pi, pi, N+1);
+      mpfr_set_z(y, pi, rnd_mode);
+      if (mpfr_cmp(x, y) != 0)
+        {
+          fprintf(stderr, "does not converge\n");
+          exit(1);
+        }
+      MPFR_SET_EXP (x, MPFR_GET_EXP(x) - 4*N);
+      mpz_clear(pi);
+      mpz_clear(num);
+      mpz_clear(den);
+      mpz_clear(d3);
+      mpz_clear(d2);
+      mpz_clear(tmp);
+      mpfr_clear(y);
     }
-    MPFR_EXP(x) -= 4*N;
-    mpz_clear(pi); mpz_clear(num); mpz_clear(den); mpz_clear(d3); mpz_clear(d2);
-    mpz_clear(tmp); mpfr_clear(y);
-  } else
-     mpfr_pi_machin3(x, rnd_mode);
+  else
+    mpfr_pi_machin3(x, rnd_mode);
   /* store computed value */
-  if (__gmpfr_const_pi_prec==0) mpfr_init2(__mpfr_const_pi, prec);
-  else mpfr_set_prec(__mpfr_const_pi, prec);
+  if (__gmpfr_const_pi_prec==0)
+    mpfr_init2(__mpfr_const_pi, prec);
+  else
+    mpfr_set_prec(__mpfr_const_pi, prec);
   mpfr_set(__mpfr_const_pi, x, rnd_mode);
   __gmpfr_const_pi_prec=prec;
   __mpfr_const_pi_rnd=rnd_mode;