Clean up code.

Add generic ZivLoop controller.


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

1 files changed, 102 insertions, 119 deletions

diff --git a/zeta.c b/zeta.c
index 958c8f8a2..18d92570d 100644
--- a/zeta.c
+++ b/zeta.c
@@ -1,6 +1,6 @@
 /* mpfr_zeta -- compute the Riemann Zeta function
 
-Copyright 2003, 2004 Free Software Foundation.
+Copyright 2003, 2004, 2005 Free Software Foundation.
 Contributed by Jean-Luc Re'my and the Spaces project, INRIA Lorraine.
 
 This file is part of the MPFR Library.
@@ -22,9 +22,7 @@ MA 02111-1307, USA. */
 
 /* #define DEBUG */
 
-#include <stdlib.h>
-#include <stdio.h>
-
+#define MPFR_NEED_LONGLONG_H
 #include "mpfr-impl.h"
 
 /*
@@ -76,11 +74,7 @@ mpfr_zeta_part_b (mpfr_t b, mpfr_srcptr s, int n, int p, mpfr_t *tc)
   mpfr_neg (s1, s1, GMP_RNDN);
   mpfr_ui_pow (u, n, s1, GMP_RNDN);
   mpfr_mul (b, d, u, GMP_RNDN);
-#ifdef DEBUG
-  printf ("\npart b=");
-  mpfr_out_str (stdout, 10, 0, b, GMP_RNDN);
-  printf ("\n");
-#endif
+  MPFR_TRACE (MPFR_DUMP (b));
   mpfr_clear (s1);
   mpfr_clear (d);
   mpfr_clear (u);
@@ -126,7 +120,7 @@ mpfr_zeta_part_a (mpfr_t sum, mpfr_srcptr s, int n)
   int i, preca;
   mpfr_t u, s1;
 
-  preca = mpfr_get_prec (sum);
+  preca = MPFR_PREC (sum);
   mpfr_init2 (u, preca);
   mpfr_init2 (s1, preca);
   mpfr_neg (s1, s, GMP_RNDN);
@@ -139,11 +133,7 @@ mpfr_zeta_part_a (mpfr_t sum, mpfr_srcptr s, int n)
       mpfr_add (sum, sum, u, GMP_RNDN);
     }
   mpfr_add_ui (sum, sum, 1, GMP_RNDN);
-#ifdef DEBUG
-  printf ("\npart a = ");
-  mpfr_out_str (stdout, 10, 0, sum, GMP_RNDN);
-  printf ("\n");
-#endif
+  MPFR_TRACE (MPFR_DUMP (sum));
   mpfr_clear (s1);
   mpfr_clear (u);
 }
@@ -156,15 +146,22 @@ mpfr_zeta_part_a (mpfr_t sum, mpfr_srcptr s, int n)
 static int
 mpfr_zeta_pos (mpfr_t z, mpfr_srcptr s, mp_rnd_t rnd_mode)
 {
-  int p, n, l, add, can_round;
+  int p, n, l, add;
   double beta, sd, dnep;
   mpfr_t a, b, c, z_pre, f, g, s1;
   mpfr_t *tc1;
   mp_prec_t precz, precs, d, dint;
   int inex;
+  MPFR_ZIV_DECL (loop);
+
+  MPFR_TRACE (MPFR_DUMP (s));
 
-  precz = mpfr_get_prec (z);
-  precs = mpfr_get_prec (s);
+  precz = MPFR_PREC (z);
+  precs = MPFR_PREC (s);
+
+  d = precz + 11;
+  /* we want that s1 = s-1 is exact, i.e. we should have PREC(s1) >= EXP(s) */
+  mpfr_init2 (s1, MAX (precs, ((mpfr_uexp_t) MPFR_EXP (s))));
 
   mpfr_init2 (a, MPFR_PREC_MIN);
   mpfr_init2 (b, MPFR_PREC_MIN);
@@ -172,22 +169,15 @@ mpfr_zeta_pos (mpfr_t z, mpfr_srcptr s, mp_rnd_t rnd_mode)
   mpfr_init2 (z_pre, MPFR_PREC_MIN);
   mpfr_init2 (f, MPFR_PREC_MIN);
   mpfr_init2 (g, MPFR_PREC_MIN);
-#ifdef DEBUG
-  printf ("Warning: mpfr_zeta assumes that s and Zeta(s) are not both representable in mpfr\n");
-  printf ("s=");
-  mpfr_print_binary (s);
-  printf ("\n");
-#endif
-  d = precz + 11;
-  /* we want that s1 = s-1 is exact, i.e. we should have PREC(s1) >= EXP(s) */
-  mpfr_init2 (s1, (precs > (mp_exp_unsigned_t) MPFR_EXP(s)) ?
-	      precs : (mp_exp_unsigned_t) MPFR_EXP(s));
-  do
+
+  MPFR_ZIV_INIT (loop, d);
+  for (;;)
     {
       /* Principal loop: we compute, in z_pre,
 	 an approximation of Zeta(s), that we send to mpfr_can_round */
       mpfr_sub_ui (s1, s, 1, GMP_RNDN);
       MPFR_ASSERTN (MPFR_IS_FP (s1));
+
       if (MPFR_IS_ZERO (s1))
         {
           mpfr_set_inf (z, 1);
@@ -199,11 +189,8 @@ mpfr_zeta_pos (mpfr_t z, mpfr_srcptr s, mp_rnd_t rnd_mode)
 	  uses the approximation Zeta(s)=1/(s-1)+gamma,
 	  where gamma is Euler's constant */
 	{
-	  dint = MAX(d + 3, precs);
-#ifdef DEBUG
-	  printf ("branch 1\n");
-	  printf ("internal precision=%d\n", dint);
-#endif
+	  dint = MAX (d + 3, precs);
+	  MPFR_TRACE (printf ("branch 1\ninternal precision=%d\n", dint));
 	  mpfr_set_prec (z_pre, dint);
 	  mpfr_set_prec (g, dint);
 	  mpfr_ui_div (z_pre, 1, s1, GMP_RNDN);
@@ -213,9 +200,8 @@ mpfr_zeta_pos (mpfr_t z, mpfr_srcptr s, mp_rnd_t rnd_mode)
       else /* Branch 2 */
         {
           size_t size;
-#ifdef DEBUG
-          printf ("branch 2\n");
-#endif
+
+          MPFR_TRACE (printf ("branch 2\n"));
           /* Computation of parameters n, p and working precision */
           dnep = (double) d * LOG2;
           sd = mpfr_get_d (s, GMP_RNDN);
@@ -235,20 +221,18 @@ mpfr_zeta_pos (mpfr_t z, mpfr_srcptr s, mp_rnd_t rnd_mode)
               p = 1 + (int) beta / 2;
               n = 1 + (int) ((sd + 2.0 * (double) p - 1.0) / 6.2832);
             }
-#ifdef DEBUG
-          printf ("\nn=%d\np=%d\n",n,p);
-#endif
+          MPFR_TRACE (printf ("\nn=%d\np=%d\n",n,p));
           /* add = 4 + floor(1.5 * log(d) / log (2)).
              We should have add >= 10, which is always fulfilled since
              d = precz + 11 >= 12, thus ceil(log2(d)) >= 4 */
-          add = 4 + (3 * __gmpfr_ceil_log2 ((double) d)) / 2;
+          add = 4 + (3 * MPFR_INT_CEIL_LOG2 (d)) / 2;
           MPFR_ASSERTD(add >= 10);
           dint = d + add;
           if (dint < precs)
             dint = precs;
-#ifdef DEBUG
-          printf("internal precision=%d\n",dint);
-#endif
+
+          MPFR_TRACE (printf("internal precision=%d\n",dint));
+
           size = (p + 1) * sizeof(mpfr_t);
           tc1 = (mpfr_t*) (*__gmp_allocate_func) (size);
           for (l=1; l<=p; l++)
@@ -258,9 +242,9 @@ mpfr_zeta_pos (mpfr_t z, mpfr_srcptr s, mp_rnd_t rnd_mode)
           mpfr_set_prec (c, dint);
           mpfr_set_prec (z_pre, dint);
           mpfr_set_prec (f, dint);
-#ifdef DEBUG
-          printf ("precision of z =%d\n", precz);
-#endif
+
+          MPFR_TRACE (printf ("precision of z =%d\n", precz));
+
           /* Computation of the coefficients c_k */
           mpfr_zeta_c (p, tc1);
           /* Computation of the 3 parts of the fonction Zeta. */
@@ -271,11 +255,7 @@ mpfr_zeta_pos (mpfr_t z, mpfr_srcptr s, mp_rnd_t rnd_mode)
           mpfr_ui_div (c, 1, s1, GMP_RNDN);
           mpfr_ui_pow (f, n, s1, GMP_RNDN);
           mpfr_div (c, c, f, GMP_RNDN);
-#ifdef DEBUG
-          printf ("\npart c=");
-          mpfr_out_str (stdout, 10, 0, c, GMP_RNDN);
-          printf ("\n");
-#endif
+	  MPFR_TRACE (MPFR_DUMP (c));
           mpfr_add (z_pre, a, c, GMP_RNDN);
           mpfr_add (z_pre, z_pre, b, GMP_RNDN);
           for (l=1; l<=p; l++)
@@ -283,34 +263,18 @@ mpfr_zeta_pos (mpfr_t z, mpfr_srcptr s, mp_rnd_t rnd_mode)
           (*__gmp_free_func) (tc1, size);
           /* End branch 2 */
         }
-#ifdef DEBUG
-      printf ("\nZeta(s) before rounding=");
-      mpfr_print_binary (z_pre);
-#endif
-      can_round = mpfr_can_round (z_pre, d - 3, GMP_RNDN, GMP_RNDZ,
-                                  precz + (rnd_mode == GMP_RNDN));
-      if (can_round)
-	{
-#ifdef DEBUG
-	  printf ("\nwe can round");
-#endif
-	}
-      else
-	{
-#ifdef DEBUG
-	  printf ("\nwe cannot round");
-#endif
-	  d = d + 5;
-	}
+
+      MPFR_TRACE (MPFR_DUMP (z_pre));
+      if (MPFR_LIKELY (mpfr_can_round (z_pre, d - 3, GMP_RNDN, GMP_RNDZ,
+				       precz + (rnd_mode == GMP_RNDN))))
+	break;
+      MPFR_ZIV_NEXT (loop, d);
     }
-  while (can_round == 0);
+  MPFR_ZIV_FREE (loop);
 
   inex = mpfr_set (z, z_pre, rnd_mode);
-#ifdef DEBUG
-  printf ("\nZeta(s) after rounding=");
-  mpfr_print_binary (z);
-  printf ("\n");
-#endif
+  MPFR_TRACE (MPFR_DUMP (z));
+
  clear_and_return:
   mpfr_clear (a);
   mpfr_clear (b);
@@ -327,81 +291,86 @@ int
 mpfr_zeta (mpfr_t z, mpfr_srcptr s, mp_rnd_t rnd_mode)
 {
   double sd, eps, m1, c;
-  int can_round;
   long add;
   mpfr_t z_pre, s1, s2, y, p;
   mp_prec_t precz, prec1, precs, precs1;
   int inex;
+  MPFR_SAVE_EXPO_DECL (expo);
 
   /* Zero, Nan or Inf ? */
-  if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(s) ))
+  if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (s)))
     {
-      if (MPFR_IS_NAN(s))
+      if (MPFR_IS_NAN (s))
 	{
 	  MPFR_SET_NAN (z);
 	  MPFR_RET_NAN;
 	}
-      else if (MPFR_IS_INF(s))
+      else if (MPFR_IS_INF (s))
 	{
-	  if (MPFR_SIGN(s) > 0)
+	  if (MPFR_IS_POS (s))
 	    return mpfr_set_ui (z, 1, GMP_RNDN); /* Zeta(+Inf) = 1 */
 	  MPFR_SET_NAN (z); /* Zeta(-Inf) = NaN */
 	  MPFR_RET_NAN;
 	}
       else /* s iz zero */
 	{
-          MPFR_ASSERTD(MPFR_IS_ZERO(s));
+          MPFR_ASSERTD (MPFR_IS_ZERO (s));
 	  mpfr_set_ui (z, 1, rnd_mode);
 	  mpfr_div_2ui (z, z, 1, rnd_mode);
-	  MPFR_CHANGE_SIGN(z);
-	  MPFR_RET(0);
+	  MPFR_CHANGE_SIGN (z);
+	  MPFR_RET (0);
 	}
     }
   MPFR_CLEAR_FLAGS(z);
 
   /* s is neither Nan, nor Inf, nor Zero */
-  mpfr_init2 (s2, mpfr_get_prec (s));
+  mpfr_init2 (s2, MPFR_PREC (s));
   mpfr_div_2ui (s2, s, 1, rnd_mode);
-  if (MPFR_IS_NEG(s) && mpfr_floor(s2, s2) == 0) /* Case s = -2n */
+  if (MPFR_IS_NEG (s) && mpfr_floor (s2, s2) == 0) /* Case s = -2n */
     {
       mpfr_clear (s2);
       return mpfr_set_ui (z, 0, rnd_mode);
     }
   mpfr_clear (s2);
-  precz = mpfr_get_prec (z);
-  precs = mpfr_get_prec (s);
-  /* Precision precs1 needed to represent 1 - s, and s + 2,
-     without any truncation */
-  precs1 = precs + 2 + MAX (0, - MPFR_GET_EXP (s));
-  sd = mpfr_get_d (s, GMP_RNDN) - 1.0;
-  if (sd < 0.0)
-    sd = -sd; /* now sd = abs(s-1.0) */
-  /* Precision prec1 is the precision on elementary computations; it ensures
-     a final precision prec1 - add for zeta(s) */
-  /* eps = pow (2.0, - (double) precz - 14.0); */
-  eps = __gmpfr_ceil_exp2 (- (double) precz - 14.0);
-  m1 = 1.0 + MAX(1.0 / eps,  2.0 * sd) * (1.0 + eps);
-  c = (1.0 + eps) * (1.0 + eps * MAX(8.0, m1));
-  /* add = 1 + floor(log(c*c*c*(13 + m1))/log(2)); */
-  add = __gmpfr_ceil_log2 (c * c * c * (13.0 + m1));
-  prec1 = precz + add; /* Note that prec1 is still incremented by  10 at the first  entry of loop below */
-  prec1 = MAX(prec1, precs1);
+
+  MPFR_SAVE_EXPO_MARK (expo);
+
+  /* Compute Zeta */
   if (MPFR_IS_POS (s) && MPFR_GET_EXP (s) >= 0) /* Case s >= 1/2 */
     inex = mpfr_zeta_pos (z, s, rnd_mode);
   else /* use reflection formula
           zeta(s) = 2^s*Pi^(s-1)*sin(Pi*s/2)*gamma(1-s)*zeta(1-s) */
     {
-      mpfr_init2 (z_pre, MPFR_PREC_MIN);
-      mpfr_init2 (s1, MPFR_PREC_MIN);
-      mpfr_init2 (y, MPFR_PREC_MIN);
-      mpfr_init2 (p, MPFR_PREC_MIN);
-      do
+      MPFR_ZIV_DECL (loop);
+
+      precz = MPFR_PREC (z);
+      precs = MPFR_PREC (s);
+
+      /* Precision precs1 needed to represent 1 - s, and s + 2,
+	 without any truncation */
+      precs1 = precs + 2 + MAX (0, - MPFR_GET_EXP (s));
+      sd = mpfr_get_d (s, GMP_RNDN) - 1.0;
+      if (sd < 0.0)
+	sd = -sd; /* now sd = abs(s-1.0) */
+      /* Precision prec1 is the precision on elementary computations;
+	 it ensures a final precision prec1 - add for zeta(s) */
+      /* eps = pow (2.0, - (double) precz - 14.0); */
+      eps = __gmpfr_ceil_exp2 (- (double) precz - 14.0);
+      m1 = 1.0 + MAX(1.0 / eps,  2.0 * sd) * (1.0 + eps);
+      c = (1.0 + eps) * (1.0 + eps * MAX(8.0, m1));
+      /* add = 1 + floor(log(c*c*c*(13 + m1))/log(2)); */
+      add = __gmpfr_ceil_log2 (c * c * c * (13.0 + m1));
+      prec1 = precz + add; 
+      prec1 = MAX (prec1, precs1) + 10;
+      
+      mpfr_init2 (z_pre, prec1);
+      mpfr_init2 (s1, prec1);
+      mpfr_init2 (y, prec1);
+      mpfr_init2 (p, prec1);
+      
+      MPFR_ZIV_INIT (loop, prec1);
+      for (;;)
 	{
-	  prec1 = prec1 + 10;
-          mpfr_set_prec (z_pre, prec1);
-          mpfr_set_prec (s1, prec1);
-          mpfr_set_prec (y, prec1);
-          mpfr_set_prec (p, prec1);
           mpfr_ui_sub (s1, 1, s, GMP_RNDN); /* s1 = 1-s */
           mpfr_zeta_pos (z_pre, s1, GMP_RNDN); /* zeta(1-s)  */
           mpfr_gamma (y, s1, GMP_RNDN);        /* gamma(1-s) */
@@ -416,15 +385,29 @@ mpfr_zeta (mpfr_t z, mpfr_srcptr s, mp_rnd_t rnd_mode)
           mpfr_pow (y, y, s1, GMP_RNDN); /* (2*Pi)^(s-1) */
           mpfr_mul (z_pre, z_pre, y, GMP_RNDN);
           mpfr_mul_2ui (z_pre, z_pre, 1, GMP_RNDN);
-          can_round = mpfr_can_round (z_pre, prec1 - add, GMP_RNDN, GMP_RNDZ,
-                                      precz + (rnd_mode == GMP_RNDN));
+
+          if (MPFR_LIKELY (mpfr_can_round (z_pre, prec1 - add, GMP_RNDN, 
+					   GMP_RNDZ,
+					   precz + (rnd_mode == GMP_RNDN))))
+	    break;
+
+	  /* Actualisation of the precision */
+	  MPFR_ZIV_NEXT (loop, prec1);
+	  mpfr_set_prec (z_pre, prec1);
+	  mpfr_set_prec (s1, prec1);
+	  mpfr_set_prec (y, prec1);
+	  mpfr_set_prec (p, prec1);
         }
-      while (can_round == 0);
+      MPFR_ZIV_FREE (loop);
+
       inex = mpfr_set (z, z_pre, rnd_mode);
+
       mpfr_clear(z_pre);
       mpfr_clear(s1);
       mpfr_clear(y);
       mpfr_clear(p);
     }
-  return inex;
+
+  MPFR_SAVE_EXPO_FREE (expo);
+  return mpfr_check_range (z, inex, rnd_mode);
 }