Clean up code.

Add generic ZivLoop controller.


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

12 files changed, 451 insertions, 421 deletions

diff --git a/gmp_op.c b/gmp_op.c
index 7ebf8a273..961071db0 100644
--- a/gmp_op.c
+++ b/gmp_op.c
@@ -143,9 +143,10 @@ int
 mpfr_add_q (mpfr_ptr y, mpfr_srcptr x, mpq_srcptr z, mp_rnd_t rnd_mode)
 {
   mpfr_t    t,q;
-  mp_prec_t p = MPFR_PREC(y)+10;
+  mp_prec_t p;
   mp_exp_t  err;
   int res;
+  MPFR_ZIV_DECL (loop);
 
   if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
     {
@@ -171,8 +172,11 @@ mpfr_add_q (mpfr_ptr y, mpfr_srcptr x, mpq_srcptr z, mp_rnd_t rnd_mode)
 	}
     }
 
+  p = MPFR_PREC (y) + 10;
   mpfr_init2 (t, p);
   mpfr_init2 (q, p);
+  
+  MPFR_ZIV_INIT (loop, p);
   for (;;)
     {
       res = mpfr_set_q (q, z, GMP_RNDN);  /* Error <= 1/2 ulp(q) */
@@ -192,18 +196,18 @@ mpfr_add_q (mpfr_ptr y, mpfr_srcptr x, mpq_srcptr z, mp_rnd_t rnd_mode)
       if (MPFR_LIKELY (!MPFR_IS_ZERO (t)))
 	{
 	  err = (mp_exp_t) p - 1 - MAX (MPFR_GET_EXP(q)-MPFR_GET_EXP(t), 0);
-	  res = mpfr_can_round (t, err, GMP_RNDN, GMP_RNDZ, 
-				MPFR_PREC(y) + (rnd_mode == GMP_RNDN) );
-	  if (MPFR_LIKELY (res != 0))  /* We can round! */
+	  if (MPFR_LIKELY (mpfr_can_round (t, err, GMP_RNDN, GMP_RNDZ, 
+				MPFR_PREC(y) + (rnd_mode == GMP_RNDN)))) 
 	    {
-	      res = mpfr_set(y, t, rnd_mode);
+	      res = mpfr_set (y, t, rnd_mode);	      
 	      break;
 	    }
 	}
-      p += BITS_PER_MP_LIMB;     /* Next precision if we continue */
+      MPFR_ZIV_NEXT (loop, p);
       mpfr_set_prec (t, p);
       mpfr_set_prec (q, p);
     } 
+  MPFR_ZIV_FREE (loop);
   mpfr_clear (t);
   mpfr_clear (q);
   return res;
@@ -213,9 +217,10 @@ int
 mpfr_sub_q (mpfr_ptr y, mpfr_srcptr x, mpq_srcptr z,mp_rnd_t rnd_mode)
 {
   mpfr_t t,q;
-  mp_prec_t p = MPFR_PREC (y)+10;
+  mp_prec_t p;
   int res;
   mp_exp_t err;
+  MPFR_ZIV_DECL (loop);
 
   if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
     {
@@ -246,8 +251,11 @@ mpfr_sub_q (mpfr_ptr y, mpfr_srcptr x, mpq_srcptr z,mp_rnd_t rnd_mode)
         }
     }
 
+  p = MPFR_PREC (y) + 10;
   mpfr_init2 (t, p);
   mpfr_init2 (q, p);
+  
+  MPFR_ZIV_INIT (loop, p);
   for(;;)
     {
       res = mpfr_set_q(q, z, GMP_RNDN);  /* Error <= 1/2 ulp(q) */
@@ -275,10 +283,11 @@ mpfr_sub_q (mpfr_ptr y, mpfr_srcptr x, mpq_srcptr z,mp_rnd_t rnd_mode)
 	      break;
 	    }
 	}
-      p += BITS_PER_MP_LIMB;             /* Next precision if we continue */
+      MPFR_ZIV_NEXT (loop, p);
       mpfr_set_prec (t, p);
       mpfr_set_prec (q, p);
     }
+  MPFR_ZIV_FREE (loop);
   mpfr_clear (t);
   mpfr_clear (q);
   return res;
diff --git a/log10.c b/log10.c
index 798b7b36f..571edc7c0 100644
--- a/log10.c
+++ b/log10.c
@@ -1,6 +1,6 @@
 /* mpfr_log10 -- logarithm in base 10.
 
-Copyright 2001, 2002, 2003, 2004 Free Software Foundation, Inc.
+Copyright 2001, 2002, 2003, 2004, 2005 Free Software Foundation, Inc.
 
 This file is part of the MPFR Library.
 
@@ -30,89 +30,81 @@ MA 02111-1307, USA. */
 int
 mpfr_log10 (mpfr_ptr r, mpfr_srcptr a, mp_rnd_t rnd_mode)
 {
-  int inexact = 0;
+  int inexact;
+  MPFR_SAVE_EXPO_DECL (expo);
 
   /* If a is NaN, the result is NaN */
-  if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(a) ))
+  if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (a)))
     {
-      if (MPFR_IS_NAN(a))
+      if (MPFR_IS_NAN (a))
 	{
-	  MPFR_SET_NAN(r);
+	  MPFR_SET_NAN (r);
 	  MPFR_RET_NAN;
 	}
       /* check for infinity before zero */
-      else if (MPFR_IS_INF(a))
+      else if (MPFR_IS_INF (a))
 	{
-	  if (MPFR_IS_NEG(a))
+	  if (MPFR_IS_NEG (a))
 	    /* log10(-Inf) = NaN */
 	    {
-	      MPFR_SET_NAN(r);
+	      MPFR_SET_NAN (r);
 	      MPFR_RET_NAN;
 	    }
 	  else /* log10(+Inf) = +Inf */
 	    {
-	      MPFR_SET_INF(r);
-	      MPFR_SET_POS(r);
-	      MPFR_RET(0); /* exact */
+	      MPFR_SET_INF (r);
+	      MPFR_SET_POS (r);
+	      MPFR_RET (0); /* exact */
 	    }
 	}
       else /* a = 0 */
 	{
-          MPFR_ASSERTD(MPFR_IS_ZERO(a));
-	  MPFR_SET_INF(r);
-	  MPFR_SET_NEG(r);
-	  MPFR_RET(0); /* log10(0) is an exact -infinity */
+          MPFR_ASSERTD (MPFR_IS_ZERO (a));
+	  MPFR_SET_INF (r);
+	  MPFR_SET_NEG (r);
+	  MPFR_RET (0); /* log10(0) is an exact -infinity */
 	}
     }
 
   /* If a is negative, the result is NaN */
-  if (MPFR_UNLIKELY( MPFR_IS_NEG(a) ))
+  if (MPFR_UNLIKELY (MPFR_IS_NEG (a)))
     {
-      MPFR_SET_NAN(r);
+      MPFR_SET_NAN (r);
       MPFR_RET_NAN;
     }
 
   /* If a is 1, the result is 0 */
   if (mpfr_cmp_ui (a, 1) == 0)
     {
-      MPFR_SET_ZERO(r);
-      MPFR_SET_POS(r);
-      MPFR_RET(0); /* result is exact */
+      MPFR_SET_ZERO (r);
+      MPFR_SET_POS (r);
+      MPFR_RET (0); /* result is exact */
     }
-
-  /* Useless due to mpfr_set
-     MPFR_CLEAR_FLAGS(r);*/
+  MPFR_SAVE_EXPO_MARK (expo);
 
   /* General case */
   {
     /* Declaration of the intermediary variable */
     mpfr_t t, tt;
-    int ok;
-
+    MPFR_ZIV_DECL (loop);
     /* Declaration of the size variable */
-    mp_prec_t Nx = MPFR_PREC(a);   /* Precision of input variable */
+    mp_prec_t Nx = MPFR_PREC(a);    /* Precision of input variable */
     mp_prec_t Ny = MPFR_PREC(r);   /* Precision of output variable */
-
-    mp_prec_t Nt;   /* Precision of the intermediary variable */
-    long int err;  /* Precision of error */
+    mp_prec_t Nt;        /* Precision of the intermediary variable */
+    mp_exp_t  err;                           /* Precision of error */
 
     /* compute the precision of intermediary variable */
-    Nt = MAX(Nx, Ny);
+    Nt = MAX (Nx, Ny);
     /* the optimal number of bits : see algorithms.ps */
     Nt = Nt + 4+ MPFR_INT_CEIL_LOG2 (Nt);
 
     /* initialise of intermediary variables */
-    mpfr_init (t);
-    mpfr_init (tt);
-
+    mpfr_init2 (t, Nt);
+    mpfr_init2 (tt, Nt);
 
     /* First computation of log10 */
-    do {
-
-      /* reactualisation of the precision */
-      mpfr_set_prec (t, Nt);
-      mpfr_set_prec (tt, Nt);
-
+    MPFR_ZIV_INIT (loop, Nt);
+    for(;;) {
       /* compute log10 */
       mpfr_set_ui (t, 10, GMP_RNDN);   /* 10 */
       mpfr_log (t, t, GMP_RNDD);       /* log(10) */
@@ -121,20 +113,25 @@ mpfr_log10 (mpfr_ptr r, mpfr_srcptr a, mp_rnd_t rnd_mode)
 
       /* estimation of the error */
       err = Nt - 4;
-
-      ok = mpfr_can_round (t, err, GMP_RNDN, GMP_RNDZ,
-                           Ny + (rnd_mode == GMP_RNDN));
-
-      /* log10(10^n) is exact */
-      if ((MPFR_SIGN(t) > 0) && mpfr_integer_p (t))
-        if (mpfr_ui_pow_ui (tt, 10, (unsigned long) (mpfr_get_d1 (t) + 0.5),
-                            GMP_RNDN) == 0)
-          if (mpfr_cmp (a, tt) == 0)
-            ok = 1;
+      if (MPFR_LIKELY (mpfr_can_round (t, err, GMP_RNDN, GMP_RNDZ,
+				       Ny + (rnd_mode == GMP_RNDN))))
+	break;
+
+      /* log10(10^n) is exact:
+	 FIXME: Can we have 10^n exactly representable as a mpfr_t
+	 but n can't fit an unsigned long? */
+      if (MPFR_IS_POS (t) 
+	  && mpfr_integer_p (t) && mpfr_fits_ulong_p (t, GMP_RNDN)
+	  && !mpfr_ui_pow_ui (tt, 10, mpfr_get_ui (t, GMP_RNDN), GMP_RNDN)
+          && mpfr_cmp (a, tt) == 0)
+	break;
 
       /* actualisation of the precision */
-      Nt += 10;
-    } while ((err < 0) || !ok);
+      MPFR_ZIV_NEXT (loop, Nt);
+      mpfr_set_prec (t, Nt);
+      mpfr_set_prec (tt, Nt);
+    }
+    MPFR_ZIV_FREE (loop);
 
     inexact = mpfr_set (r, t, rnd_mode);
 
@@ -142,5 +139,6 @@ mpfr_log10 (mpfr_ptr r, mpfr_srcptr a, mp_rnd_t rnd_mode)
     mpfr_clear (tt);
   }
 
-  return inexact;
+  MPFR_SAVE_EXPO_FREE (expo);
+  return mpfr_check_range (r, inexact, rnd_mode);
 }
diff --git a/log1p.c b/log1p.c
index 400f9b590..fb3d298f6 100644
--- a/log1p.c
+++ b/log1p.c
@@ -1,6 +1,6 @@
 /* mpfr_log1p -- Compute log(1+x)
 
-Copyright 2001, 2002, 2003, 2004 Free Software Foundation, Inc.
+Copyright 2001, 2002, 2003, 2004, 2005 Free Software Foundation, Inc.
 
 This file is part of the MPFR Library.
 
@@ -28,36 +28,37 @@ MA 02111-1307, USA. */
 int
 mpfr_log1p (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode)
 {
-  int comp, inexact = 0;
+  int comp, inexact;
+  MPFR_SAVE_EXPO_DECL (expo);
 
-  if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(x)))
+  if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
     {
-      if (MPFR_IS_NAN(x)) 
+      if (MPFR_IS_NAN (x)) 
 	{
-	  MPFR_SET_NAN(y); 
+	  MPFR_SET_NAN (y); 
 	  MPFR_RET_NAN;
 	}
       /* check for inf or -inf (result is not defined) */
-      else if (MPFR_IS_INF(x))
+      else if (MPFR_IS_INF (x))
 	{
-	  if (MPFR_IS_POS(x))
+	  if (MPFR_IS_POS (x))
 	    {
-	      MPFR_SET_INF(y);
-	      MPFR_SET_POS(y);
-	      MPFR_RET(0);
+	      MPFR_SET_INF (y);
+	      MPFR_SET_POS (y);
+	      MPFR_RET (0);
 	    }
 	  else
 	    {
-	      MPFR_SET_NAN(y);
+	      MPFR_SET_NAN (y);
 	      MPFR_RET_NAN;
 	    }
 	}
       else /* x is zero */
 	{
-          MPFR_ASSERTD(MPFR_IS_ZERO(x));
-	  MPFR_SET_ZERO(y);   /* log1p(+/- 0) = +/- 0 */
-	  MPFR_SET_SAME_SIGN(y, x);
-	  MPFR_RET(0);
+          MPFR_ASSERTD (MPFR_IS_ZERO (x));
+	  MPFR_SET_ZERO (y);   /* log1p(+/- 0) = +/- 0 */
+	  MPFR_SET_SAME_SIGN (y, x);
+	  MPFR_RET (0);
 	}
     }
   
@@ -68,60 +69,61 @@ mpfr_log1p (mpfr_ptr y, mpfr_srcptr x, mp_rnd_t rnd_mode)
       if (comp == 0)
 	/* x=0: log1p(-1)=-inf (division by zero) */
 	{
-	  MPFR_SET_INF(y);
-	  MPFR_SET_NEG(y);
-	  MPFR_RET(0);
+	  MPFR_SET_INF (y);
+	  MPFR_SET_NEG (y);
+	  MPFR_RET (0);
 	}
-      MPFR_SET_NAN(y);
+      MPFR_SET_NAN (y);
       MPFR_RET_NAN;
     }
 
-  MPFR_CLEAR_FLAGS(y);
+  MPFR_SAVE_EXPO_MARK (expo);
 
   /* General case */
   {
     /* Declaration of the intermediary variable */
     mpfr_t t;
-
     /* Declaration of the size variable */
     mp_prec_t Nx = MPFR_PREC(x);   /* Precision of input variable */
-    mp_prec_t Ny = MPFR_PREC(y);   /* Precision of input variable */
-    
-    mp_prec_t Nt;   /* Precision of the intermediary variable */
-    long int err;  /* Precision of error */
-                
+    mp_prec_t Ny = MPFR_PREC(y);   /* Precision of input variable */    
+    mp_prec_t Nt;       /* Precision of the intermediary variable */
+    mp_exp_t err;                           /* Precision of error */
+    MPFR_ZIV_DECL (loop);
+
     /* compute the precision of intermediary variable */
-    Nt = MAX(Nx,Ny);
+    Nt = MAX (Nx, Ny);
     /* the optimal number of bits : see algorithms.ps */
     Nt = Nt + 5 + MPFR_INT_CEIL_LOG2 (Nt);
 
     /* initialise of intermediary variable */
-    mpfr_init (t);
+    mpfr_init2 (t, Nt);
 
     /* First computation of cosh */
-    do
+    MPFR_ZIV_INIT (loop, Nt);
+    for (;;)
       {
-        /* reactualisation of the precision */
-        mpfr_set_prec (t, Nt);
-      
         /* compute log1p */
-        mpfr_add_ui (t, x, 1, GMP_RNDN);   /* 1+x */
+        mpfr_add_ui (t, x, 1, GMP_RNDN);      /* 1+x */
         mpfr_log (t, t, GMP_RNDN);        /* log(1+x)*/
 
         /* estimation of the error */
         /*err=Nt-(__gmpfr_ceil_log2(1+pow(2,1-MPFR_GET_EXP(t))));*/
         err = Nt - (MAX (1 - MPFR_GET_EXP (t), 0) + 1);
 
+	if (MPFR_LIKELY (mpfr_can_round (t, err, GMP_RNDN, GMP_RNDZ,
+                                         Ny + (rnd_mode == GMP_RNDN))))
+	  break;
+	
         /* actualisation of the precision */
-        Nt += 10;
+	MPFR_ZIV_NEXT (loop, Nt);
+	mpfr_set_prec (t, Nt);
       }
-    while ((err < 0) || !mpfr_can_round (t, err, GMP_RNDN, GMP_RNDZ,
-                                         Ny + (rnd_mode == GMP_RNDN)));
-
+    MPFR_ZIV_FREE (loop);
     inexact = mpfr_set (y, t, rnd_mode);
 
     mpfr_clear (t);
   }
 
-  return inexact;
+  MPFR_SAVE_EXPO_FREE (expo);
+  return mpfr_check_range (y, inexact, rnd_mode);
 }
diff --git a/log2.c b/log2.c
index 548cb123d..173acd92f 100644
--- a/log2.c
+++ b/log2.c
@@ -1,6 +1,6 @@
 /* mpfr_log2 -- log base 2
 
-Copyright 2001, 2002, 2003, 2004 Free Software Foundation, Inc.
+Copyright 2001, 2002, 2003, 2004, 2005 Free Software Foundation, Inc.
 
 This file is part of the MPFR Library.
 
@@ -23,110 +23,109 @@ MA 02111-1307, USA. */
 #include "mpfr-impl.h"
 
  /* The computation of r=log2(a)
-
-    r=log2(a)=log(a)/log(2)
- */
+      r=log2(a)=log(a)/log(2)      */
 
 int
 mpfr_log2 (mpfr_ptr r, mpfr_srcptr a, mp_rnd_t rnd_mode)
 {
-  int inexact = 0;
+  int inexact;
+  MPFR_SAVE_EXPO_DECL (expo);
 
-  if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(a) ))
+  if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (a)))
     {
       /* If a is NaN, the result is NaN */
-      if (MPFR_IS_NAN(a))
+      if (MPFR_IS_NAN (a))
 	{
-	  MPFR_SET_NAN(r);
+	  MPFR_SET_NAN (r);
 	  MPFR_RET_NAN;
 	}
       /* check for infinity before zero */
-      else if (MPFR_IS_INF(a))
+      else if (MPFR_IS_INF (a))
 	{
-	  if (MPFR_IS_NEG(a))
+	  if (MPFR_IS_NEG (a))
 	    /* log(-Inf) = NaN */
 	    {
-	      MPFR_SET_NAN(r);
+	      MPFR_SET_NAN (r);
 	      MPFR_RET_NAN;
 	    }
 	  else /* log(+Inf) = +Inf */
 	    {
-	      MPFR_SET_INF(r);
-	      MPFR_SET_POS(r);
-	      MPFR_RET(0);
+	      MPFR_SET_INF (r);
+	      MPFR_SET_POS (r);
+	      MPFR_RET (0);
 	    }
 	}
       else /* a is zero */
 	{
-          MPFR_ASSERTD(MPFR_IS_ZERO(a));
-	  MPFR_SET_INF(r);
-	  MPFR_SET_NEG(r);
-	  MPFR_RET(0); /* log2(0) is an exact -infinity */
+          MPFR_ASSERTD (MPFR_IS_ZERO (a));
+	  MPFR_SET_INF (r);
+	  MPFR_SET_NEG (r);
+	  MPFR_RET (0); /* log2(0) is an exact -infinity */
 	}
     }
 
   /* If a is negative, the result is NaN */
-  if (MPFR_UNLIKELY(MPFR_IS_NEG(a)))
+  if (MPFR_UNLIKELY (MPFR_IS_NEG (a)))
     {
-      MPFR_SET_NAN(r);
+      MPFR_SET_NAN (r);
       MPFR_RET_NAN;
     }
 
   /* If a is 1, the result is 0 */
-  if (mpfr_cmp_ui(a, 1) == 0)
+  if (MPFR_UNLIKELY (mpfr_cmp_ui (a, 1) == 0))
     {
-      MPFR_SET_ZERO(r);
-      MPFR_SET_POS(r);
-      MPFR_RET(0); /* only "normal" case where the result is exact */
+      MPFR_SET_ZERO (r);
+      MPFR_SET_POS (r);
+      MPFR_RET (0); /* only "normal" case where the result is exact */
     }
 
-  /* If a is integer, log2(a) is exact*/
-  if (mpfr_cmp_ui_2exp (a, 1, MPFR_GET_EXP (a) - 1) == 0)
+  /* If a is 2^N, log2(a) is exact*/
+  if (MPFR_UNLIKELY (mpfr_cmp_ui_2exp (a, 1, MPFR_GET_EXP (a) - 1) == 0))
     return mpfr_set_si(r, MPFR_GET_EXP (a) - 1, rnd_mode);
 
+  MPFR_SAVE_EXPO_MARK (expo);
+
   /* General case */
   {
     /* Declaration of the intermediary variable */
     mpfr_t t, tt;
-
     /* Declaration of the size variable */
     mp_prec_t Nx = MPFR_PREC(a);   /* Precision of input variable */
     mp_prec_t Ny = MPFR_PREC(r);   /* Precision of input variable */
-
-    mp_prec_t Nt;   /* Precision of the intermediary variable */
-    long int err;  /* Precision of error */
-
+    mp_prec_t Nt;       /* Precision of the intermediary variable */
+    mp_exp_t err;                           /* Precision of error */
+    MPFR_ZIV_DECL (loop);
 
     /* compute the precision of intermediary variable */
-    Nt=MAX(Nx,Ny);
+    Nt = MAX (Nx, Ny);
     /* the optimal number of bits : see algorithms.ps */
-    Nt=Nt + 3 + MPFR_INT_CEIL_LOG2 (Nt);
+    Nt = Nt + 3 + MPFR_INT_CEIL_LOG2 (Nt);
 
     /* initialise of intermediary	variable */
-    mpfr_init(t);
-    mpfr_init(tt);
-
+    mpfr_init2 (t, Nt);
+    mpfr_init2 (tt, Nt);
 
     /* First computation of log2 */
-    do
+    MPFR_ZIV_INIT (loop, Nt);
+    for (;;)
       {
-        /* reactualisation of the precision */
-        mpfr_set_prec(t,Nt);
-        mpfr_set_prec(tt,Nt);
-
         /* compute log2 */
         mpfr_const_log2(t,GMP_RNDD); /* log(2) */
         mpfr_log(tt,a,GMP_RNDN);     /* log(a) */
         mpfr_div(t,tt,t,GMP_RNDN); /* log(a)/log(2) */
-
+	
         /* estimation of the error */
-        err=Nt-3;
-
+        err = Nt-3;
+	if (MPFR_LIKELY (mpfr_can_round (t, err, GMP_RNDN, GMP_RNDZ,
+                                         Ny + (rnd_mode == GMP_RNDN))))
+	  break;
+	
         /* actualisation of the precision */
-        Nt += 10;
+	MPFR_ZIV_NEXT (loop, Nt);
+        mpfr_set_prec (t, Nt);
+        mpfr_set_prec (tt, Nt);
       }
-    while ((err < 0) || !mpfr_can_round (t, err, GMP_RNDN, GMP_RNDZ,
-                                         Ny + (rnd_mode == GMP_RNDN)));
+    MPFR_ZIV_FREE (loop);
 
     inexact = mpfr_set (r, t, rnd_mode);
 
@@ -134,5 +133,6 @@ mpfr_log2 (mpfr_ptr r, mpfr_srcptr a, mp_rnd_t rnd_mode)
     mpfr_clear (tt);
   }
 
-  return inexact;
+  MPFR_SAVE_EXPO_FREE (expo);
+  return mpfr_check_range (r, inexact, rnd_mode);
 }
diff --git a/pow.c b/pow.c
index a750d2f9a..2729ea2a3 100644
--- a/pow.c
+++ b/pow.c
@@ -199,8 +199,9 @@ int
 mpfr_pow (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y, mp_rnd_t rnd_mode)
 {
   int inexact = 1;
+  MPFR_SAVE_EXPO_DECL (expo);
 
-  if (MPFR_ARE_SINGULAR(x,y))
+  if (MPFR_ARE_SINGULAR (x, y))
     {
       /* pow(x, 0) returns 1 for any x, even a NaN. */
       if (MPFR_UNLIKELY( MPFR_IS_ZERO(y) ))
@@ -336,69 +337,63 @@ mpfr_pow (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y, mp_rnd_t rnd_mode)
       MPFR_RET_NAN;
     }
 
-  MPFR_CLEAR_FLAGS(z);
+  MPFR_SAVE_EXPO_MARK (expo);
 
   /* General case */
   {
     /* Declaration of the intermediary variable */
-    mpfr_t t, te, ti;
-    int loop = 0, ok;
-
+    mpfr_t t;
+    int loop = 0;
     /* Declaration of the size variable */
     mp_prec_t Nx = MPFR_PREC(x);   /* Precision of input variable */
     mp_prec_t Ny = MPFR_PREC(y);   /* Precision of input variable */
     mp_prec_t Nz = MPFR_PREC(z);   /* Precision of output variable */
-
-    mp_prec_t Nt;   /* Precision of the intermediary variable */
-    long int err;  /* Precision of error */
+    mp_prec_t Nt;        /* Precision of the intermediary variable */
+    mp_exp_t err, exp_te;                    /* Precision of error */
+    MPFR_ZIV_DECL (ziv_loop);
 
     /* compute the precision of intermediary variable */
-    Nt = MAX(Nx, Ny);
-    Nt = MAX(Nt, Nz); /* take account of the output precision too! */
+    Nt = MAX (Nx, Ny);
+    Nt = MAX (Nt, Nz); /* take account of the output precision too! */
     /* the optimal number of bits : see algorithms.ps */
     Nt = Nt + 5 + MPFR_INT_CEIL_LOG2 (Nt);
 
     /* initialise of intermediary variable */
-    mpfr_init (t);
-    mpfr_init (ti);
-    mpfr_init (te);
+    mpfr_init2 (t, Nt);
 
-    do
+    MPFR_ZIV_INIT (ziv_loop, Nt);
+    for (;;)
       {
         loop ++;
 
-        /* reactualisation of the precision */
-        mpfr_set_prec (t, Nt);
-        mpfr_set_prec (ti, Nt);
-        mpfr_set_prec (te, Nt);
-
         /* compute exp(y*ln(x)) */
-        mpfr_log (ti, x, GMP_RNDU);         /* ln(x) */
-        mpfr_mul (te, y, ti, GMP_RNDU);       /* y*ln(x) */
-        mpfr_exp (t, te, GMP_RNDN);         /* exp(y*ln(x))*/
-
-        MPFR_ASSERTN(MPFR_IS_FP(te));
-        MPFR_ASSERTN(MPFR_NOTZERO(te));
-        /* otherwise MPFR_EXP(te) below doesn't exist */
+        mpfr_log (t, x, GMP_RNDU);               /* ln(x) */
+        mpfr_mul (t, y, t, GMP_RNDU);            /* y*ln(x) */
+	exp_te = MPFR_GET_EXP (t);               /* FIXME: May overflow */
+        mpfr_exp (t, t, GMP_RNDN);               /* exp(y*ln(x))*/
+                                                 /* FIXME: May overflow */
 
 	/* estimate of the error -- see pow function in algorithms.ps */
-        err = Nt - (MPFR_GET_EXP (te) + 3);
-
-	/* actualisation of the precision */
-        Nt += 10;
-
-        ok = mpfr_can_round (t, err, GMP_RNDN, GMP_RNDZ,
-                             Nz + (rnd_mode == GMP_RNDN));
+        err = Nt - (exp_te + 3);
+        if (MPFR_LIKELY (mpfr_can_round (t, err, GMP_RNDN, GMP_RNDZ,
+					 Nz + (rnd_mode == GMP_RNDN))))
+	  break;
 
         /* check exact power */
-        if (ok == 0 && loop == 1)
+        if (loop == 1)
           {
-            ok = mpfr_pow_is_exact (x, y);
-            if (ok)
-              inexact = 0;
-          }
+            if (mpfr_pow_is_exact (x, y))
+	      {
+		inexact = 0;
+		break;
+	      }
+	  }
+
+        /* reactualisation of the precision */
+	MPFR_ZIV_NEXT (ziv_loop, Nt);
+        mpfr_set_prec (t, Nt);
       }
-    while (err < 0 || ok == 0);
+    MPFR_ZIV_FREE (ziv_loop);
 
     if (inexact)
       inexact = mpfr_set (z, t, rnd_mode);
@@ -406,9 +401,8 @@ mpfr_pow (mpfr_ptr z, mpfr_srcptr x, mpfr_srcptr y, mp_rnd_t rnd_mode)
       mpfr_set (z, t, GMP_RNDN);
 
     mpfr_clear (t);
-    mpfr_clear (ti);
-    mpfr_clear (te);
   }
 
-  return inexact;
+  MPFR_SAVE_EXPO_FREE (expo);
+  return mpfr_check_range (z, inexact, rnd_mode);
 }
diff --git a/pow_si.c b/pow_si.c
index 8e2ccfa2a..fc4cd3f26 100644
--- a/pow_si.c
+++ b/pow_si.c
@@ -35,7 +35,7 @@ mpfr_pow_si (mpfr_ptr y, mpfr_srcptr x, long int n, mp_rnd_t rnd_mode)
     return mpfr_pow_ui (y, x, n, rnd_mode);
   else
     {
-      if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(x) ))
+      if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
 	{
 	  if (MPFR_IS_NAN (x))
 	    {
@@ -53,7 +53,7 @@ mpfr_pow_si (mpfr_ptr y, mpfr_srcptr x, long int n, mp_rnd_t rnd_mode)
 	    }
 	  else /* x is zero */
 	    {
-              MPFR_ASSERTD (MPFR_IS_ZERO(x));
+              MPFR_ASSERTD (MPFR_IS_ZERO (x));
 	      MPFR_SET_INF(y);
 	      if (MPFR_IS_POS (x) || ((unsigned) n & 1) == 0)
 		MPFR_SET_POS (y);
@@ -62,7 +62,7 @@ mpfr_pow_si (mpfr_ptr y, mpfr_srcptr x, long int n, mp_rnd_t rnd_mode)
 	      MPFR_RET(0);
 	    }
 	}
-      MPFR_CLEAR_FLAGS(y);
+      MPFR_CLEAR_FLAGS (y);
 
       /* detect exact powers: x^(-n) is exact iff x is a power of 2 */
       if (mpfr_cmp_si_2exp (x, MPFR_SIGN(x), MPFR_EXP(x) - 1) == 0)
@@ -73,13 +73,13 @@ mpfr_pow_si (mpfr_ptr y, mpfr_srcptr x, long int n, mp_rnd_t rnd_mode)
 	  expx --;
 	  MPFR_ASSERTD (n < 0);
 	  /* Warning n*expx may overflow! 
-	     Many systems abort with LONG_MIN / 1 or LONG_MIN/-1*/
+	     Some systems abort with LONG_MIN / 1 or LONG_MIN/-1*/
 	  if (n != -1 && expx > 0 && -expx < MPFR_EXP_MIN / (-n))
 	    MPFR_EXP (y) = MPFR_EMIN_MIN - 1; /* Underflow */
 	  else if (n != -1 && expx < 0 && -expx > MPFR_EXP_MAX / (-n))
 	    MPFR_EXP (y) = MPFR_EMAX_MAX + 1; /* Overflow */
           else
-	    MPFR_EXP(y) += n * expx;
+	    MPFR_EXP (y) += n * expx;
           return mpfr_check_range (y, 0, rnd_mode);
         }
 
@@ -89,18 +89,17 @@ mpfr_pow_si (mpfr_ptr y, mpfr_srcptr x, long int n, mp_rnd_t rnd_mode)
       {
         /* Declaration of the intermediary variable */
         mpfr_t t;
-
         /* Declaration of the size variable */
-        mp_prec_t Nx = MPFR_PREC(x);   /* Precision of input variable */
-        mp_prec_t Ny = MPFR_PREC(y);   /* Precision of output variable */
-
-        mp_prec_t Nt;   /* Precision of the intermediary variable */
-        long int err;   /* Precision of error */
+        mp_prec_t Nx = MPFR_PREC (x);   /* Precision of input variable */
+        mp_prec_t Ny = MPFR_PREC (y);   /* Precision of output variable */
+        mp_prec_t Nt;         /* Precision of the intermediary variable */
+        mp_exp_t  err;                            /* Precision of error */
         int inexact;
 	MPFR_SAVE_EXPO_DECL (expo);
+	MPFR_ZIV_DECL (loop);
 
         /* compute the precision of intermediary variable */
-        Nt = MAX(Nx,Ny);
+        Nt = MAX (Nx, Ny);
         /* the optimal number of bits : see algorithms.ps */
         Nt = Nt + 3 + MPFR_INT_CEIL_LOG2 (Nt);
 
@@ -108,24 +107,28 @@ mpfr_pow_si (mpfr_ptr y, mpfr_srcptr x, long int n, mp_rnd_t rnd_mode)
 
         /* initialise of intermediary	variable */
         mpfr_init2 (t, Nt);
-
-        do {
-	  /* reactualisation of the precision */
-	  mpfr_set_prec (t, Nt);
-	  
-	  /* compute 1/(x^n) n>0*/
-	  mpfr_pow_ui (t, x, (unsigned long int) n, GMP_RNDN);
-	  inexact = MPFR_IS_ZERO (t) || MPFR_IS_INF (t);
-	  mpfr_ui_div (t, 1, t, GMP_RNDN);
-	  inexact = inexact || MPFR_IS_ZERO (t) || MPFR_IS_INF (t);
-	  /* error estimate -- see pow function in algorithms.ps */
-	  err = Nt - 3;
-	  
-	  /* actualisation of the precision */
-	  Nt += BITS_PER_MP_LIMB;
-	} while (inexact == 0 &&
-		 !mpfr_can_round (t, err, GMP_RNDN, GMP_RNDZ,
-				  Ny + (rnd_mode == GMP_RNDN)));
+	
+	MPFR_ZIV_INIT (loop, Nt);
+        for (;;)
+	  {
+	    /* compute 1/(x^n) n>0*/
+	    mpfr_pow_ui (t, x, (unsigned long int) n, GMP_RNDN);
+	    inexact = MPFR_IS_ZERO (t) || MPFR_IS_INF (t);
+	    mpfr_ui_div (t, 1, t, GMP_RNDN);
+	    inexact = inexact || MPFR_IS_ZERO (t) || MPFR_IS_INF (t);
+
+	    /* error estimate -- see pow function in algorithms.ps */
+	    err = Nt - 3;
+	    if (MPFR_LIKELY (inexact != 0
+			     || mpfr_can_round (t, err, GMP_RNDN, GMP_RNDZ,
+						Ny + (rnd_mode == GMP_RNDN))))
+	      break;
+
+	    /* actualisation of the precision */
+	    Nt += BITS_PER_MP_LIMB;
+	    mpfr_set_prec (t, Nt);
+	  }
+	MPFR_ZIV_FREE (loop);
 
         inexact = mpfr_set (y, t, rnd_mode);
         mpfr_clear (t);
diff --git a/pow_ui.c b/pow_ui.c
index 17d3a97fd..097968c83 100644
--- a/pow_ui.c
+++ b/pow_ui.c
@@ -33,6 +33,7 @@ mpfr_pow_ui (mpfr_ptr x, mpfr_srcptr y, unsigned long int n, mp_rnd_t rnd)
   int inexact;
   mp_rnd_t rnd1;
   MPFR_SAVE_EXPO_DECL (expo);
+  MPFR_ZIV_DECL (loop);
 
   if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (y)))
     {
@@ -81,41 +82,48 @@ mpfr_pow_ui (mpfr_ptr x, mpfr_srcptr y, unsigned long int n, mp_rnd_t rnd)
 	/* y^2 = sqr(y) */
 	return mpfr_sqr (x, y, rnd);
     }
-  /* MPFR_CLEAR_FLAGS useless due to mpfr_set */
 
+  /* Augment exponent range */
   MPFR_SAVE_EXPO_MARK (expo);
   __gmpfr_emin -= 3;  /* So that we can check for underflow properly */
 
-  prec = MPFR_PREC (x);
-  mpfr_init2 (res, prec + 9);
+  /* setup initial precision */
+  prec = MPFR_PREC (x) + 3 + BITS_PER_MP_LIMB;
+  mpfr_init2 (res, prec);
 
   rnd1 = MPFR_IS_POS (y) ? GMP_RNDU : GMP_RNDD; /* away */
 
-  do {
-    int i;
-    prec += 3;
-    for (m = n, i = 0; m; i++, m >>= 1)
-      prec++;
-    /* now 2^(i-1) <= n < 2^i */
-    mpfr_set_prec (res, prec);
-    err = prec <= (mpfr_prec_t) i ? 0 : prec - (mpfr_prec_t) i;
-    MPFR_ASSERTD (i >= 1);
-    mpfr_clear_overflow ();
-    mpfr_clear_underflow ();
-    /* First step: compute square from y */
-    inexact = mpfr_sqr (res, y, GMP_RNDU);
-    if (n & (1UL << (i-2)))
-      inexact |= mpfr_mul (res, res, y, rnd1);
-    for (i -= 3; i >= 0 && !mpfr_overflow_p () && !mpfr_underflow_p (); i--) 
-      {
-	inexact |= mpfr_sqr (res, res, GMP_RNDU);
-	if (n & (1UL << i))
-	  inexact |= mpfr_mul (res, res, y, rnd1);
-      }
-  }
-  while (inexact && !mpfr_overflow_p () && !mpfr_underflow_p ()
-	 && !mpfr_can_round (res, err, GMP_RNDN, GMP_RNDZ,
-			     MPFR_PREC(x) + (rnd == GMP_RNDN)));
+  MPFR_ZIV_INIT (loop, prec);
+  for (;;) 
+    {
+      int i;
+      for (m = n, i = 0; m; i++, m >>= 1);
+      /* now 2^(i-1) <= n < 2^i */
+      err = prec <= (mpfr_prec_t) i ? 0 : prec - (mpfr_prec_t) i;
+      MPFR_ASSERTD (i >= 1);
+      mpfr_clear_overflow ();
+      mpfr_clear_underflow ();
+      /* First step: compute square from y */
+      inexact = mpfr_sqr (res, y, GMP_RNDU);
+      if (n & (1UL << (i-2)))
+	inexact |= mpfr_mul (res, res, y, rnd1);
+      for (i -= 3; i >= 0 && !mpfr_overflow_p () && !mpfr_underflow_p (); i--) 
+	{
+	  inexact |= mpfr_sqr (res, res, GMP_RNDU);
+	  if (n & (1UL << i))
+	    inexact |= mpfr_mul (res, res, y, rnd1);
+	}
+      if (MPFR_LIKELY (inexact == 0 
+		       || mpfr_overflow_p () || mpfr_underflow_p ()
+		       || mpfr_can_round (res, err, GMP_RNDN, GMP_RNDZ,
+					  MPFR_PREC(x) + (rnd == GMP_RNDN))))
+	break;
+      /* Actualisation of the precision */
+      MPFR_ZIV_NEXT (loop, prec);
+      mpfr_set_prec (res, prec);      
+    }
+  MPFR_ZIV_FREE (loop);
+
   inexact = mpfr_set (x, res, rnd);
   mpfr_clear (res);
 
diff --git a/pow_z.c b/pow_z.c
index 1f324e4a7..14c0be8cf 100644
--- a/pow_z.c
+++ b/pow_z.c
@@ -30,31 +30,32 @@ mpfr_pow_pos_z (mpfr_ptr y, mpfr_srcptr x, mpz_srcptr z, mp_rnd_t rnd)
   int inexact;
   mp_rnd_t rnd1;
   mpz_t absz;
+  mp_size_t size_z;
+  MPFR_ZIV_DECL (loop);
 
   MPFR_ASSERTD (mpz_sgn (z) != 0);
 
   if (MPFR_UNLIKELY (mpz_cmpabs_ui (z, 1) == 0))
     return mpfr_set (y, x, rnd);
 
-  prec = MPFR_PREC (x);
-  mpfr_init2 (res, prec + 9);
-
   rnd1 = MPFR_IS_POS (x) ? GMP_RNDU : GMP_RNDD; /* away */
   absz[0] = z[0];
   SIZ (absz) = ABS(SIZ(absz)); /* Hack to get abs(z) */
+  MPFR_MPZ_SIZEINBASE2 (size_z, z);
 
-  do {
-    mp_size_t i;
-    
-    MPFR_MPZ_SIZEINBASE2 (i, z);
+  prec = MPFR_PREC (x) + 3 + size_z;
+  mpfr_init2 (res, prec);
+
+  MPFR_ZIV_INIT (loop, prec);
+  for (;;) {
+    mp_size_t i = size_z;
     /* now 2^(i-1) <= z < 2^i */
-    prec += 3 + i;
-    mpfr_set_prec (res, prec);
     err = prec <= (mpfr_prec_t) i ? 0 : prec - (mpfr_prec_t) i;
     MPFR_ASSERTD (i >= 2);
     mpfr_clear_overflow ();
     mpfr_clear_underflow ();
-    /* First step: compute square from y */
+ 
+   /* First step: compute square from y */
     inexact = mpfr_mul (res, x, x, GMP_RNDU);
     if (mpz_tstbit (absz, i-2))
       inexact |= mpfr_mul (res, res, x, rnd1);
@@ -63,11 +64,20 @@ mpfr_pow_pos_z (mpfr_ptr y, mpfr_srcptr x, mpz_srcptr z, mp_rnd_t rnd)
 	inexact |= mpfr_sqr (res, res, GMP_RNDU);
 	if (mpz_tstbit (absz, i))
 	  inexact |= mpfr_mul (res, res, x, rnd1);      
-      } /* for */
-  } while (inexact && !mpfr_overflow_p() && !mpfr_underflow_p ()
-	   && !mpfr_can_round (res, err, GMP_RNDN, GMP_RNDZ,
-			       MPFR_PREC (y) + (rnd == GMP_RNDN)));
-  
+      }
+
+    if (MPFR_LIKELY (inexact == 0 
+		     || mpfr_overflow_p () || mpfr_underflow_p ()
+		     || mpfr_can_round (res, err, GMP_RNDN, GMP_RNDZ,
+					MPFR_PREC (y) + (rnd == GMP_RNDN))))
+      break;
+
+    /* Actualisation of the precision */
+    MPFR_ZIV_NEXT (loop, prec);
+    mpfr_set_prec (res, prec);
+  }
+  MPFR_ZIV_FREE (loop);
+
   inexact = mpfr_set (y, res, rnd);
   mpfr_clear (res);
 
@@ -189,7 +199,8 @@ mpfr_pow_z (mpfr_ptr y, mpfr_srcptr x, mpz_srcptr z, mp_rnd_t rnd)
       /* Declaration of the intermediary variable */
       mpfr_t t;
       mp_prec_t Nt;   /* Precision of the intermediary variable */
-      
+      MPFR_ZIV_DECL (loop);
+
       /* compute the precision of intermediary variable */
       Nt = MAX (MPFR_PREC (x), MPFR_PREC (y));
       
@@ -198,7 +209,8 @@ mpfr_pow_z (mpfr_ptr y, mpfr_srcptr x, mpz_srcptr z, mp_rnd_t rnd)
     
       /* initialise of intermediary	variable */
       mpfr_init2 (t, Nt);
-      
+
+      MPFR_ZIV_INIT (loop, Nt);
       for (;;) {
 	/* compute 1/(x^n) n>0 */
 	mpfr_pow_pos_z (t, x, z, GMP_RNDN);
@@ -206,15 +218,16 @@ mpfr_pow_z (mpfr_ptr y, mpfr_srcptr x, mpz_srcptr z, mp_rnd_t rnd)
 	mpfr_ui_div (t, 1, t, GMP_RNDN);
 	inexact = inexact || MPFR_IS_ZERO (t) || MPFR_IS_INF (t);
 
-	if (inexact != 0
-	    || mpfr_can_round (t, Nt - 3, GMP_RNDN, GMP_RNDZ,
-			       MPFR_PREC (y) + (rnd == GMP_RNDN)))
+	if (MPFR_LIKELY (inexact != 0
+			 || mpfr_can_round (t, Nt - 3, GMP_RNDN, GMP_RNDZ,
+					 MPFR_PREC (y) + (rnd == GMP_RNDN))))
 	  break;
 	/* actualisation of the precision */
-	Nt += BITS_PER_MP_LIMB;
+	MPFR_ZIV_NEXT (loop, Nt);
 	mpfr_set_prec (t, Nt);
       }
-      
+      MPFR_ZIV_FREE (loop);
+
       inexact = mpfr_set (y, t, rnd);
       mpfr_clear (t);
     }
diff --git a/strtofr.c b/strtofr.c
index 70c191652..5a01b5712 100644
--- a/strtofr.c
+++ b/strtofr.c
@@ -409,6 +409,7 @@ parsed_string_to_mpfr (mpfr_t x, struct parsed_string *pstr, mp_rnd_t rnd)
   size_t pstr_size;
   mp_size_t ysize, real_ysize;
   int res, err;
+  MPFR_ZIV_DECL (loop);
   TMP_DECL (marker);
 
   /* determine the minimal precision for the computation */
@@ -416,6 +417,7 @@ parsed_string_to_mpfr (mpfr_t x, struct parsed_string *pstr, mp_rnd_t rnd)
 
   /* compute y as long as rounding is not possible */
   TMP_MARK(marker);
+  MPFR_ZIV_INIT (loop, prec);
   for (;;)
     {
       /* initialize y to the value of 0.mant_s[0]...mant_s[pr-1] */
@@ -633,8 +635,9 @@ parsed_string_to_mpfr (mpfr_t x, struct parsed_string *pstr, mp_rnd_t rnd)
         break;
 
       /* update the prec for next loop */
-      prec += BITS_PER_MP_LIMB;
+      MPFR_ZIV_NEXT (loop, prec);
     } /* loop */
+  MPFR_ZIV_FREE (loop);
 
   /* round y */
   if (mpfr_round_raw (MPFR_MANT (x), result,
diff --git a/sum.c b/sum.c
index ecde90397..650a7feca 100644
--- a/sum.c
+++ b/sum.c
@@ -193,13 +193,13 @@ static int mpfr_list_sum_once (mpfr_ptr ret, mpfr_srcptr const tab[],
 int mpfr_sum (mpfr_ptr ret, mpfr_ptr const tab[], unsigned long n, 
               mp_rnd_t rnd)
 {
-  mp_prec_t initial_f, current_f;
+  mp_prec_t prec;
   unsigned long k;
   mpfr_srcptr *perm;
-  unsigned int guard_digits;
-  unsigned int initial_guard_digits;
   int error_trap;
   mpfr_t cur_sum;
+  MPFR_ZIV_DECL (loop);
+  MPFR_SAVE_EXPO_DECL (expo);
   TMP_DECL(marker);
     
   TMP_MARK(marker);
@@ -209,30 +209,39 @@ int mpfr_sum (mpfr_ptr ret, mpfr_ptr const tab[], unsigned long n,
     return 0;
   }
 
+  /* Sort */
   perm = (mpfr_srcptr *) TMP_ALLOC(n * sizeof(mpfr_srcptr)); 
-
   mpfr_count_sort (tab, n, perm);
 
-  initial_f = MAX (MPFR_PREC(tab[0]), MPFR_PREC(ret));
+  /* Initial precision */
+  prec = MAX (MPFR_PREC (tab[0]), MPFR_PREC (ret));
   k = MPFR_INT_CEIL_LOG2 (n) + 1;
-  mpfr_init2 (cur_sum, initial_f);
-  initial_guard_digits = k + 2;
-  guard_digits = initial_guard_digits;
-  do
-  {
-      current_f = initial_f + guard_digits;
-      mpfr_set_prec (cur_sum, current_f);
-      error_trap = mpfr_list_sum_once (cur_sum, perm, n, 
-                                       current_f + k);
-      guard_digits *= 2;
+  prec += k + 2;
+  mpfr_init2 (cur_sum, prec);
+
+  MPFR_SAVE_EXPO_MARK (expo);
+
+  /* Ziv Loop */
+  MPFR_ZIV_INIT (loop, prec);
+  for (;;) {
+    error_trap = mpfr_list_sum_once (cur_sum, perm, n, prec + k);
+    if (MPFR_LIKELY (error_trap == 0
+		     || mpfr_can_round (cur_sum, 
+					MPFR_GET_EXP (cur_sum) - prec + 2,
+					GMP_RNDN, rnd, MPFR_PREC (ret))))
+      break;
+    MPFR_ZIV_NEXT (loop, prec);
+    mpfr_set_prec (cur_sum, prec);
   }
-  while ((error_trap != 0) &&
-          !(mpfr_can_round (cur_sum, MPFR_GET_EXP(cur_sum) - current_f + 2,
-                            GMP_RNDN, rnd, MPFR_PREC(ret))));
+  MPFR_ZIV_FREE (loop);
+
   error_trap |= mpfr_set (ret, cur_sum, rnd);
   mpfr_clear (cur_sum);
   TMP_FREE(marker);
-  return error_trap;
+
+  MPFR_SAVE_EXPO_FREE (expo);
+  error_trap |= mpfr_check_range (ret, 0, rnd);
+  return error_trap; /* It doesn't return the ternary value */
 }
 
 
diff --git a/ui_pow_ui.c b/ui_pow_ui.c
index 478d94f2d..1303a1884 100644
--- a/ui_pow_ui.c
+++ b/ui_pow_ui.c
@@ -1,6 +1,7 @@
 /*  mpfr_ui_pow_ui -- compute the power beetween two machine integer
 
-Copyright 1999, 2000, 2001, 2002, 2003, 2004 Free Software Foundation, Inc.
+Copyright 1999, 2000, 2001, 2002, 2003, 2004, 2005 
+  Free Software Foundation, Inc.
 
 This file is part of the MPFR Library.
 
@@ -25,15 +26,15 @@ int
 mpfr_ui_pow_ui (mpfr_ptr x, unsigned long int y, unsigned long int n,
                 mp_rnd_t rnd)
 {
-  long int err;
+  mp_exp_t err;
   unsigned long m;
   mpfr_t res;
   mp_prec_t prec;
+  int size_n;
   int inexact;
+  MPFR_ZIV_DECL (loop);
   MPFR_SAVE_EXPO_DECL (expo);
 
-  MPFR_CLEAR_FLAGS(x);
-
   if (MPFR_UNLIKELY (n == 0))
     /* y^0 = 1 for any y */
     return mpfr_set_ui (x, 1, rnd);
@@ -42,18 +43,17 @@ mpfr_ui_pow_ui (mpfr_ptr x, unsigned long int y, unsigned long int n,
     /* 0^n = 0 for any n > 0 */
     return mpfr_set_ui (x, 0, rnd);
 
+  for (size_n = 0, m = n; m; size_n++, m >>= 1);
+
   MPFR_SAVE_EXPO_MARK (expo);
-  prec = MPFR_PREC (x);
-  mpfr_init2 (res, 2*prec);
+  prec = MPFR_PREC (x) + 3 + size_n;
+  mpfr_init2 (res, prec);
 
-  do
+  MPFR_ZIV_INIT (loop, prec);
+  for (;;)
     {
-      int i;
+      int i = size_n;
 
-      prec += 3;
-      for (i = 0, m = n; m; i++, m >>= 1)
-        prec++;
-      mpfr_set_prec (res, prec);
       inexact = mpfr_set_ui (res, y, GMP_RNDU);
       err = 1;
       /* now 2^(i-1) <= n < 2^i: i=1+floor(log2(n)) */
@@ -68,9 +68,17 @@ mpfr_ui_pow_ui (mpfr_ptr x, unsigned long int y, unsigned long int n,
          we have err = 1+floor(log2(n)).
          Since prec >= MPFR_PREC(x) + 4 + floor(log2(n)), prec > err */
       err = prec - err;
+      
+      if (MPFR_LIKELY (inexact == 0
+		       || mpfr_can_round (res, err, GMP_RNDN, GMP_RNDZ,
+					  MPFR_PREC(x) + (rnd == GMP_RNDN))))
+	break;
+      
+      /* Actualisation of the precision */
+      MPFR_ZIV_NEXT (loop, prec);
+      mpfr_set_prec (res, prec);
     }
-  while (inexact && !mpfr_can_round (res, err, GMP_RNDN, GMP_RNDZ,
-                                     MPFR_PREC(x) + (rnd == GMP_RNDN)));
+  MPFR_ZIV_FREE (loop);
 
   inexact = mpfr_set (x, res, rnd);
 
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);
 }