Clean up code (Fix integer types + rewrite some code).

Optimize the code by reducing the number of used variables inside the loop.


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

1 files changed, 74 insertions, 127 deletions

diff --git a/atan.c b/atan.c
index f1f87b1bc..fa43e6828 100644
--- a/atan.c
+++ b/atan.c
@@ -23,11 +23,9 @@ MA 02111-1307, USA. */
 #include "mpfr-impl.h"
 
 #define CST   2.27  /* CST=1+ln(2.4)/ln(2) */
-#define CST2  1.45  /* CST2=1/ln(2) */
 
 static int mpfr_atan_aux _MPFR_PROTO((mpfr_ptr, mpz_srcptr, long, int));
 
-#undef B
 #define A
 #define A1 1
 #define A2 2
@@ -47,213 +45,162 @@ static int mpfr_atan_aux _MPFR_PROTO((mpfr_ptr, mpz_srcptr, long, int));
 #undef GENERIC
 
 int
-mpfr_atan (mpfr_ptr arctangent, mpfr_srcptr x, mp_rnd_t rnd_mode)
+mpfr_atan (mpfr_ptr atan, mpfr_srcptr x, mp_rnd_t rnd_mode)
 {
-  mpfr_t Pisur2;
   mpfr_t xp;
-  mpfr_t arctgt;
-
-  int comparaison, sign, supplement, inexact;
-
-  mpfr_t t_arctan;
-  int i;
-  mpz_t ukz;
-  mpfr_t ukf;
-  mpfr_t sk,Ak;
-  mpz_t square;
-  mpfr_t tmp_arctan;
-  mpfr_t tmp, tmp2;
-#ifdef DEBUG
-  mpfr_t tst;
-#endif
-  int twopoweri;
-  int Prec;
-  int prec_x;
-  int prec_arctan;
-  int realprec;
-  int estimated_delta;
-  /* calculation of the floor */
+  mpfr_t arctgt, sk, tmp, tmp2;
+  mpz_t  ukz;
+  int comparaison, sign, inexact;
+  mp_exp_t supplement, estimated_delta;
+  mp_prec_t prec, realprec;
   mp_exp_t exptol;
-  int N0;
-  int logn;
+  int i, n0;
+  unsigned long twopoweri;
   MPFR_SAVE_EXPO_DECL (expo);
 
-  /* Trivial cases */
-  if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(x) ))
+  /* Singular cases */
+  if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
     {
-      if (MPFR_IS_NAN(x))
+      if (MPFR_IS_NAN (x))
 	{
-	  MPFR_SET_NAN(arctangent);
+	  MPFR_SET_NAN (atan);
 	  MPFR_RET_NAN;
 	}
-      else if (MPFR_IS_INF(x))
+      else if (MPFR_IS_INF (x))
 	{
-	  if (MPFR_IS_POS(x))
-	    /* arctan(+inf) = Pi/2 */
-	    inexact = mpfr_const_pi (arctangent, rnd_mode);
-	  else 
-	    /* arctan(-inf) = -Pi/2 */
+	  if (MPFR_IS_POS (x))  /* arctan(+inf) = Pi/2 */
+	    inexact = mpfr_const_pi (atan, rnd_mode);
+	  else /* arctan(-inf) = -Pi/2 */
 	    {
-	      inexact = -mpfr_const_pi (arctangent, 
+	      inexact = -mpfr_const_pi (atan, 
 					MPFR_INVERT_RND (rnd_mode));
-	      MPFR_CHANGE_SIGN (arctangent);
+	      MPFR_CHANGE_SIGN (atan);
 	    }
-	  MPFR_SET_EXP (arctangent, MPFR_GET_EXP (arctangent) - 1);
+	  mpfr_div_2ui (atan, atan, 1, rnd_mode);
 	  return inexact;
 	}
       else /* x is necessarily 0 */
 	{
           MPFR_ASSERTD (MPFR_IS_ZERO (x));
-	  mpfr_set_ui (arctangent, 0, GMP_RNDN);
-	  return 0; /* exact result */
+	  return mpfr_set_ui (atan, 0, rnd_mode);
 	}
     }
-  MPFR_CLEAR_FLAGS(arctangent);
-
-  sign = MPFR_SIGN(x);
-  prec_arctan = MPFR_PREC(arctangent);
 
   /* Set x_p=|x| */
-  mpfr_init2 (xp, MPFR_PREC(x));
-  mpfr_abs (xp, x, rnd_mode);
+  sign = MPFR_SIGN (x);
+  xp[0] = x[0];                 /* Hack to avoid copying X */
+  MPFR_SET_POS (xp);
 
   /* Other simple case arctang(-+1)=-+pi/4 */
   comparaison = mpfr_cmp_ui (xp, 1);
-  if (MPFR_UNLIKELY(comparaison == 0))
+  if (MPFR_UNLIKELY (comparaison == 0))
     {
-      inexact = mpfr_const_pi (arctangent, MPFR_IS_POS_SIGN(sign) ? rnd_mode
-                               : MPFR_INVERT_RND(rnd_mode));
-      MPFR_SET_EXP (arctangent, MPFR_GET_EXP (arctangent) - 2);
+      inexact = mpfr_const_pi (atan, MPFR_IS_POS_SIGN (sign) ? rnd_mode
+                               : MPFR_INVERT_RND (rnd_mode));
+      mpfr_div_2ui (atan, atan, 2, rnd_mode);
       if (MPFR_IS_NEG_SIGN (sign))
         {
           inexact = -inexact;
-          MPFR_CHANGE_SIGN (arctangent);
+          MPFR_CHANGE_SIGN (atan);
         }
-      mpfr_clear (xp);
       return inexact;
     }
 
-  if (comparaison > 0)
-    supplement = 2;
-  else
-    supplement = 2 - MPFR_GET_EXP (xp);
+  supplement = 2 - (comparaison > 0) ? 0 : MPFR_GET_EXP (xp);
 
   MPFR_SAVE_EXPO_MARK (expo);
 
   /* FIXME: Could use MPFR_INT_CEIL_LOG2? */
-  prec_x = __gmpfr_ceil_log2 ((double) MPFR_PREC(x) / BITS_PER_MP_LIMB);
-  logn   = __gmpfr_ceil_log2 (prec_x);
-  if (logn < 2) 
-    logn = 2;
-  realprec = prec_arctan + MPFR_INT_CEIL_LOG2 (prec_arctan) + 4;
+  prec = __gmpfr_ceil_log2 ((double) MPFR_PREC(x) / BITS_PER_MP_LIMB);
+
+  realprec = MPFR_PREC (atan) + MPFR_INT_CEIL_LOG2 (MPFR_PREC (atan)) + 4;
 
   mpz_init (ukz);
-  mpz_init (square);
 
   /* Initialisation    */
   mpfr_init (sk);
-  mpfr_init (ukf);
-  mpfr_init (t_arctan);
-  mpfr_init (tmp_arctan);
-  mpfr_init (tmp);
   mpfr_init (tmp2);
-  mpfr_init (Ak);
+  mpfr_init (tmp);
   mpfr_init (arctgt);
-  mpfr_init (Pisur2);
 
-  while (1)
+  for (;;)
     {
-      N0 = __gmpfr_ceil_log2 ((double) realprec + supplement + CST);
-      estimated_delta = 1 + supplement + MPFR_INT_CEIL_LOG2 (3*N0-2);
-      Prec = realprec + estimated_delta;
-
-      /* Initialisation    */
-      mpfr_set_prec (sk,Prec);
-      mpfr_set_prec (ukf, Prec);
-      mpfr_set_prec (t_arctan, Prec);
-      mpfr_set_prec (tmp_arctan, Prec);
-      mpfr_set_prec (tmp, Prec);
-      mpfr_set_prec (tmp2, Prec);
-      mpfr_set_prec (Ak, Prec);
-      mpfr_set_prec (arctgt, Prec);
+      n0 = __gmpfr_ceil_log2 ((double) realprec + supplement + CST);
+      MPFR_ASSERTD (3*n0 > 2);
+      estimated_delta = 1 + supplement + MPFR_INT_CEIL_LOG2 (3*n0-2);
+      prec = realprec + estimated_delta;
+
+      /* Initialisation */
+      mpfr_set_prec (sk, prec);
+      mpfr_set_prec (tmp2, prec);
+      mpfr_set_prec (tmp, prec);
+      mpfr_set_prec (arctgt, prec);
 
       if (comparaison > 0)
-        {
-          mpfr_set_prec (Pisur2, Prec);
-          mpfr_const_pi (Pisur2, GMP_RNDN);
-          mpfr_div_2ui (Pisur2, Pisur2, 1, GMP_RNDN);
-          mpfr_ui_div (sk, 1, xp, GMP_RNDN);
-        }
+	mpfr_ui_div (sk, 1, xp, GMP_RNDN);
       else
 	mpfr_set (sk, xp, GMP_RNDN);
 
       /* sk is 1/|x| if |x| > 1, and |x| otherwise, i.e. min(|x|, 1/|x|) */
 
       /* Assignation  */
-      mpfr_set_ui (tmp_arctan, 0, GMP_RNDN);
+      mpfr_set_ui (arctgt, 0, GMP_RNDN);
       twopoweri = 1;
-      for (i = 0; i <= N0; i++)
+      for (i = 0; i <= n0; i++)
         {
           mpfr_mul_2ui (tmp, sk, twopoweri, GMP_RNDN);
+
           /* Calculation of  trunc(tmp) --> mpz */
-          mpfr_trunc (ukf, tmp);
-          exptol = mpfr_get_z_exp (ukz, ukf);
+          mpfr_trunc (tmp, tmp);
+          exptol = mpfr_get_z_exp (ukz, tmp);
           /* since the s_k are decreasing (see algorithms.tex),
              and s_0 = min(|x|, 1/|x|) < 1, we have sk < 1,
              thus exptol < 0 */
-          MPFR_ASSERTD(exptol < 0);
+          MPFR_ASSERTD (exptol < 0);
           mpz_tdiv_q_2exp (ukz, ukz, (unsigned long int) (-exptol));
 
           /* Calculation of arctan(Ak) */
-          mpz_mul (square, ukz, ukz);
-	  mpz_neg (square, square);
-          mpfr_atan_aux (t_arctan, square, 2*twopoweri, N0 - i);
-          mpfr_set_z (Ak, ukz, GMP_RNDN);
-          mpfr_div_2ui (Ak, Ak, twopoweri, GMP_RNDN);
-          mpfr_mul (t_arctan, t_arctan, Ak, GMP_RNDN);
+	  mpfr_set_z (tmp, ukz, GMP_RNDN);
+	  mpfr_div_2ui (tmp, tmp, twopoweri, GMP_RNDN);
+          mpz_mul (ukz, ukz, ukz);
+	  mpz_neg (ukz, ukz);
+          mpfr_atan_aux (tmp2, ukz, 2*twopoweri, n0 - i);
+          mpfr_mul (tmp2, tmp2, tmp, GMP_RNDN);
 
           /* Addition and iteration */
-          mpfr_add (tmp_arctan, tmp_arctan, t_arctan, GMP_RNDN);
-          if (i < N0)
+          mpfr_add (arctgt, arctgt, tmp2, GMP_RNDN);
+          if (MPFR_LIKELY (i < n0))
             {
-              mpfr_sub (tmp, sk, Ak, GMP_RNDN);
-              mpfr_mul (tmp2, sk, Ak, GMP_RNDN);
-              mpfr_add_ui (tmp2, tmp2, 1, GMP_RNDN);
-              mpfr_div (sk, tmp, tmp2, GMP_RNDN);
+              mpfr_sub (tmp2, sk, tmp, GMP_RNDN);
+              mpfr_mul (sk, sk, tmp, GMP_RNDN);
+              mpfr_add_ui (sk, sk, 1, GMP_RNDN);
+              mpfr_div (sk, tmp2, sk, GMP_RNDN);
               twopoweri <<= 1;
             }
         }
 
       if (comparaison > 0)
-	mpfr_sub(arctgt, Pisur2, tmp_arctan, GMP_RNDN);
-      else
-	mpfr_set(arctgt, tmp_arctan, GMP_RNDN);
-      MPFR_SET_POS(arctgt);
+	{
+	  mpfr_const_pi (tmp, GMP_RNDN);
+	  mpfr_div_2ui (tmp, tmp, 1, GMP_RNDN);
+	  mpfr_sub (arctgt, tmp, arctgt, GMP_RNDN);
+	}
+      MPFR_SET_POS (arctgt);
 
-      if (!mpfr_can_round (arctgt, realprec, GMP_RNDN, GMP_RNDZ,
-                          MPFR_PREC (arctangent) + (rnd_mode == GMP_RNDN)))
-	realprec += MPFR_INT_CEIL_LOG2 (realprec);
-      else
+      if (mpfr_can_round (arctgt, realprec, GMP_RNDN, GMP_RNDZ,
+                          MPFR_PREC (atan) + (rnd_mode == GMP_RNDN)))
 	break;
+      realprec += MPFR_INT_CEIL_LOG2 (realprec);
     }
 
-  inexact = MPFR_IS_POS_SIGN(sign) ? mpfr_set (arctangent, arctgt, rnd_mode)
-    : mpfr_neg (arctangent, arctgt, rnd_mode);
+  inexact = mpfr_set4 (atan, arctgt, rnd_mode, sign);
 
   mpfr_clear (sk);
-  mpfr_clear (ukf);
-  mpfr_clear (t_arctan);
-  mpfr_clear (tmp_arctan);
-  mpfr_clear (tmp);
   mpfr_clear (tmp2);
-  mpfr_clear (Ak);
+  mpfr_clear (tmp);
   mpfr_clear (arctgt);
 
-  mpfr_clear (Pisur2);
-
-  mpfr_clear (xp);
   mpz_clear (ukz);
-  mpz_clear (square);
 
   MPFR_SAVE_EXPO_FREE (expo);
   return mpfr_check_range (arctgt, inexact, rnd_mode);