Optimize mpfr_asin by improving the choice of the initial precision.

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

1 files changed, 21 insertions, 21 deletions

diff --git a/asin.c b/asin.c
index a580cde3a..23dbebba1 100644
--- a/asin.c
+++ b/asin.c
@@ -24,8 +24,8 @@ MA 02111-1307, USA. */
 int
 mpfr_asin (mpfr_ptr asin, mpfr_srcptr x, mp_rnd_t rnd_mode)
 {
-  mpfr_t xp, arcs, tmp;
-  int sign, compared, inexact;
+  mpfr_t xp;
+  int compared, inexact;
   mp_prec_t prec;
   mp_exp_t xp_exp;
   MPFR_SAVE_EXPO_DECL (expo);
@@ -41,14 +41,12 @@ mpfr_asin (mpfr_ptr asin, mpfr_srcptr x, mp_rnd_t rnd_mode)
       else /* x = 0 */
 	{
           MPFR_ASSERTD (MPFR_IS_ZERO (x));
-	  mpfr_set_ui (asin, 0, GMP_RNDN);
+	  MPFR_SET_ZERO (asin);
 	  MPFR_RET (0); /* exact result */
 	}
     }
-  MPFR_CLEAR_FLAGS (asin);
 
   /* 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);
@@ -65,7 +63,7 @@ mpfr_asin (mpfr_ptr asin, mpfr_srcptr x, mp_rnd_t rnd_mode)
 	}
       else                              /* x = 1 or x = -1 */
 	{
-	  if (MPFR_IS_POS_SIGN (sign)) /* asin(+1) = Pi/2 */
+	  if (MPFR_IS_POS (x)) /* asin(+1) = Pi/2 */
 	    inexact = mpfr_const_pi (asin, rnd_mode);
 	  else /* asin(-1) = -Pi/2 */
 	    {
@@ -82,36 +80,38 @@ mpfr_asin (mpfr_ptr asin, mpfr_srcptr x, mp_rnd_t rnd_mode)
   /* Compute exponent of 1 - ABS(x) */
   mpfr_ui_sub (xp, 1, xp, GMP_RNDD);
   MPFR_ASSERTD (MPFR_GET_EXP (xp) <= 0);
+  MPFR_ASSERTD (MPFR_GET_EXP (x) <= 0);
   xp_exp = 2 - MPFR_GET_EXP (xp); 
-  mpfr_clear (xp);
 
   /* Set up initial prec */
   prec = MPFR_PREC (asin) + 10 + xp_exp;
-  mpfr_init2 (tmp, prec);
-  mpfr_init2 (arcs, prec);
+
+  /* If x ~ 2^-N, asin(x) ~ x + x^3/6
+     If Prec < 2*N, we can't round since x^3/6 won't be count. */
+  if (MPFR_PREC (asin) >= MPFR_PREC (x) 
+      && prec <= -2*MPFR_GET_EXP (x) + 10)
+    prec = -2*MPFR_GET_EXP (x) + 10;
 
   /* use asin(x) = atan(x/sqrt(1-x^2)) */
+
   for (;;)
     {
-      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, prec - xp_exp, GMP_RNDN, GMP_RNDZ,
+      mpfr_set_prec (xp, prec);
+      mpfr_sqr (xp, x, GMP_RNDN);
+      mpfr_ui_sub (xp, 1, xp, GMP_RNDN);
+      mpfr_sqrt (xp, xp, GMP_RNDN);
+      mpfr_div (xp, x, xp, GMP_RNDN);
+      mpfr_atan (xp, xp, GMP_RNDN);
+      if (mpfr_can_round (xp, prec - xp_exp, GMP_RNDN, GMP_RNDZ,
                           MPFR_PREC (asin) + (rnd_mode == GMP_RNDN)))
 	break;
       prec += BITS_PER_MP_LIMB;
-      mpfr_set_prec (tmp, prec);
-      mpfr_set_prec (arcs, prec);
     }
 
-  inexact = mpfr_set (asin, arcs, rnd_mode);
+  inexact = mpfr_set (asin, xp, rnd_mode);
 
-  mpfr_clear (tmp);
-  mpfr_clear (arcs);
+  mpfr_clear (xp);
 
   MPFR_SAVE_EXPO_FREE (expo);
-
   return mpfr_check_range (asin, inexact, rnd_mode);
 }