Clean up the code (Removing useless variables and avoid mixing wrong integer types).

Optimize the code by improving memory allocation scheme and by incrementating by
BITS_PER_MP_LIMB in case of an error instead of MPFR_INT_CEIL_LOG2 (prec).
Fix an overflow bug for X=+/-1 (If PI is inside the exponent range, but not PI/2).


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

1 files changed, 47 insertions, 56 deletions

diff --git a/asin.c b/asin.c
index 19c9e3975..85914ecef 100644
--- a/asin.c
+++ b/asin.c
@@ -19,106 +19,97 @@ along with the MPFR Library; see the file COPYING.LIB.  If not, write to
 the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
 MA 02111-1307, USA. */
 
-#define MPFR_NEED_LONGLONG_H
 #include "mpfr-impl.h"
 
 int
 mpfr_asin (mpfr_ptr asin, mpfr_srcptr x, mp_rnd_t rnd_mode)
 {
-  mpfr_t xp;
-  mpfr_t arcs;
-  int sign, supplement;
-  mpfr_t tmp;
-  int Prec;
-  int prec_asin;
-  int realprec;
-  int estimated_delta;
-  int compared;
-  int inexact;
+  mpfr_t xp, arcs, tmp;
+  int sign, compared, inexact;
+  mp_prec_t prec;
+  mp_exp_t xp_exp;
   MPFR_SAVE_EXPO_DECL (expo);
 
   /* Special cases */
-  if (MPFR_UNLIKELY( MPFR_IS_SINGULAR(x) ))
+  if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x)))
     {
-      if (MPFR_IS_NAN(x) || MPFR_IS_INF(x))
+      if (MPFR_IS_NAN (x) || MPFR_IS_INF (x))
 	{
-	  MPFR_SET_NAN(asin);
+	  MPFR_SET_NAN (asin);
 	  MPFR_RET_NAN;
 	}
       else /* x = 0 */
 	{
-          MPFR_ASSERTD(MPFR_IS_ZERO(x));
+          MPFR_ASSERTD (MPFR_IS_ZERO (x));
 	  mpfr_set_ui (asin, 0, GMP_RNDN);
-	  MPFR_RET(0); /* exact result */
+	  MPFR_RET (0); /* exact result */
 	}
-      MPFR_ASSERTN(0);
     }
-  MPFR_CLEAR_FLAGS(asin);
+  MPFR_CLEAR_FLAGS (asin);
 
-  /* Set x_p=|x| */
-  sign = MPFR_SIGN(x);
-  mpfr_init2 (xp, MPFR_PREC(x));
-  mpfr_abs (xp, x, rnd_mode);
+  /* Set x_p=|x| (x is a normal number) */
+  sign = MPFR_SIGN (x);
+  mpfr_init2 (xp, MPFR_PREC (x));
+  inexact = mpfr_abs (xp, x, GMP_RNDN); 
+  MPFR_ASSERTD (inexact == 0);
 
   compared = mpfr_cmp_ui (xp, 1);
 
-  if (MPFR_UNLIKELY(compared > 0))     /* asin(x) = NaN for |x| > 1 */
+  if (MPFR_UNLIKELY (compared >= 0))
     {
-      MPFR_SET_NAN(asin);
       mpfr_clear (xp);
-      MPFR_RET_NAN;
-    }
-
-  if (MPFR_UNLIKELY(compared == 0))   /* x = 1 or x = -1 */
-    {
-      if (MPFR_IS_POS_SIGN(sign)) /* asin(+1) = Pi/2 */
-        inexact = mpfr_const_pi (asin, rnd_mode);
-      else /* asin(-1) = -Pi/2 */
-        {
-          inexact = -mpfr_const_pi (asin, MPFR_INVERT_RND(rnd_mode));
-          mpfr_neg (asin, asin, rnd_mode);
-        }
-      MPFR_SET_EXP (asin, MPFR_GET_EXP (asin) - 1);
-      mpfr_clear (xp);
-      return inexact;
+      if (compared > 0)                  /* asin(x) = NaN for |x| > 1 */
+	{
+	  MPFR_SET_NAN (asin);
+	  MPFR_RET_NAN;
+	}
+      else                              /* x = 1 or x = -1 */
+	{
+	  if (MPFR_IS_POS_SIGN (sign)) /* asin(+1) = Pi/2 */
+	    inexact = mpfr_const_pi (asin, rnd_mode);
+	  else /* asin(-1) = -Pi/2 */
+	    {
+	      inexact = -mpfr_const_pi (asin, MPFR_INVERT_RND(rnd_mode));
+	      mpfr_neg (asin, asin, rnd_mode);
+	    }
+	  mpfr_div_2ui (asin, asin, 1, rnd_mode); /* May underflow */
+	  return inexact;
+	}
     }
 
   MPFR_SAVE_EXPO_MARK (expo);
 
-  prec_asin = MPFR_PREC(asin);
+  /* Compute exponent of 1 - ABS(x) */
   mpfr_ui_sub (xp, 1, xp, GMP_RNDD);
+  MPFR_ASSERTD (MPFR_GET_EXP (xp) <= 0);
+  xp_exp = 2 - MPFR_GET_EXP (xp); 
+  mpfr_clear (xp);
 
-  supplement = 2 - MPFR_GET_EXP (xp);
-  realprec = prec_asin + 10;
-
-  mpfr_init (tmp);
-  mpfr_init (arcs);
+  /* Set up initial prec */
+  prec = MPFR_PREC (asin) + 10 + xp_exp;
+  mpfr_init2 (tmp, prec);
+  mpfr_init2 (arcs, prec);
 
   /* use asin(x) = atan(x/sqrt(1-x^2)) */
   for (;;)
     {
-      estimated_delta = 1 + supplement;
-      Prec = realprec+estimated_delta;
-
-      /* Fix prec   */
-      mpfr_set_prec (tmp, Prec);
-      mpfr_set_prec (arcs, Prec);
-      mpfr_mul (tmp, x, x, GMP_RNDN);
+      mpfr_sqr (tmp, x, GMP_RNDN);
       mpfr_ui_sub (tmp, 1, tmp, GMP_RNDN);
       mpfr_sqrt (tmp, tmp, GMP_RNDN);
       mpfr_div (tmp, x, tmp, GMP_RNDN);
       mpfr_atan (arcs, tmp, GMP_RNDN);
-      if (mpfr_can_round (arcs, realprec, GMP_RNDN, GMP_RNDZ,
-                          MPFR_PREC(asin) + (rnd_mode == GMP_RNDN)))
+      if (mpfr_can_round (arcs, prec - xp_exp, GMP_RNDN, GMP_RNDZ,
+                          MPFR_PREC (asin) + (rnd_mode == GMP_RNDN)))
 	break;
-      realprec += MPFR_INT_CEIL_LOG2 (realprec);
+      prec += BITS_PER_MP_LIMB;
+      mpfr_set_prec (tmp, prec);
+      mpfr_set_prec (arcs, prec);
     }
 
   inexact = mpfr_set (asin, arcs, rnd_mode);
 
   mpfr_clear (tmp);
   mpfr_clear (arcs);
-  mpfr_clear (xp);
 
   MPFR_SAVE_EXPO_FREE (expo);