tests/tatan.c: added testcase for bug found by Christopher Creutzig

               (atan2_different_prec). Check underflow of arctan(Inf) and
               arctan(1) with reduced exponent range.
atan.c: fixed bug when x is very near but differs from 1 atan.c: expand
        exponent range when computing arctan(Inf) and arctan(+-1).
[merge r6180 r6181 r6186 r6187 from trunk]


git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/branches/2.4@6214 280ebfd0-de03-0410-8827-d642c229c3f4

1 files changed, 27 insertions, 12 deletions

diff --git a/atan.c b/atan.c
index 21144ba54..e9e64ec18 100644
--- a/atan.c
+++ b/atan.c
@@ -192,7 +192,7 @@ mpfr_atan (mpfr_ptr atan, mpfr_srcptr x, mp_rnd_t rnd_mode)
   mp_exp_t exptol;
   mp_prec_t prec, realprec;
   unsigned long twopoweri;
-  int comparaison, inexact, inexact2;
+  int comparaison, inexact;
   int i, n0, oldn0;
   MPFR_GROUP_DECL (group);
   MPFR_SAVE_EXPO_DECL (expo);
@@ -211,6 +211,7 @@ mpfr_atan (mpfr_ptr atan, mpfr_srcptr x, mp_rnd_t rnd_mode)
         }
       else if (MPFR_IS_INF (x))
         {
+          MPFR_SAVE_EXPO_MARK (expo);
           if (MPFR_IS_POS (x))  /* arctan(+inf) = Pi/2 */
             inexact = mpfr_const_pi (atan, rnd_mode);
           else /* arctan(-inf) = -Pi/2 */
@@ -219,10 +220,9 @@ mpfr_atan (mpfr_ptr atan, mpfr_srcptr x, mp_rnd_t rnd_mode)
                                         MPFR_INVERT_RND (rnd_mode));
               MPFR_CHANGE_SIGN (atan);
             }
-          inexact2 = mpfr_div_2ui (atan, atan, 1, rnd_mode);
-          if (MPFR_UNLIKELY (inexact2))
-            inexact = inexact2; /* An underflow occurs */
-          MPFR_RET (inexact);
+          mpfr_div_2ui (atan, atan, 1, rnd_mode);  /* exact (no exceptions) */
+          MPFR_SAVE_EXPO_FREE (expo);
+          return mpfr_check_range (atan, inexact, rnd_mode);
         }
       else /* x is necessarily 0 */
         {
@@ -242,6 +242,8 @@ mpfr_atan (mpfr_ptr atan, mpfr_srcptr x, mp_rnd_t rnd_mode)
   /* Set x_p=|x| */
   MPFR_TMP_INIT_ABS (xp, x);
 
+  MPFR_SAVE_EXPO_MARK (expo);
+
   /* Other simple case arctan(-+1)=-+pi/4 */
   comparaison = mpfr_cmp_ui (xp, 1);
   if (MPFR_UNLIKELY (comparaison == 0))
@@ -254,17 +256,14 @@ mpfr_atan (mpfr_ptr atan, mpfr_srcptr x, mp_rnd_t rnd_mode)
           inexact = -inexact;
           MPFR_CHANGE_SIGN (atan);
         }
-      inexact2 = mpfr_div_2ui (atan, atan, 2, rnd_mode);
-      if (MPFR_UNLIKELY (inexact2))
-        inexact = inexact2; /* an underflow occurs */
-      return inexact;
+      mpfr_div_2ui (atan, atan, 2, rnd_mode);  /* exact (no exceptions) */
+      MPFR_SAVE_EXPO_FREE (expo);
+      return mpfr_check_range (atan, inexact, rnd_mode);
     }
 
   realprec = MPFR_PREC (atan) + MPFR_INT_CEIL_LOG2 (MPFR_PREC (atan)) + 4;
   prec = realprec + BITS_PER_MP_LIMB;
 
-  MPFR_SAVE_EXPO_MARK (expo);
-
   /* Initialisation    */
   mpz_init (ukz);
   MPFR_GROUP_INIT_4 (group, prec, sk, tmp, tmp2, arctgt);
@@ -314,9 +313,24 @@ mpfr_atan (mpfr_ptr atan, mpfr_srcptr x, mp_rnd_t rnd_mode)
 
       /* sk is 1/|x| if |x| > 1, and |x| otherwise, i.e. min(|x|, 1/|x|) */
 
+      /* If sk=1, then if |x| < 1, we have 1 - 2^(-prec-1) <= |x| < 1,
+         or if |x| > 1, we have 1 - 2^(-prec-1) <= 1/|x| < 1, thus in all
+         cases ||x| - 1| <= 2^(-prec), from which it follows
+         |atan|x| - Pi/4| <= 2^(-prec), given the Taylor expansion
+         atan(1+x) = Pi/4 + x/2 - x^2/4 + ...
+         Since Pi/4 = 0.785..., the error is at most one ulp.
+      */
+      if (MPFR_UNLIKELY(mpfr_cmp_ui (sk, 1) == 0))
+        {
+          mpfr_const_pi (arctgt, GMP_RNDN); /* 1/2 ulp extra error */
+          mpfr_div_2ui (arctgt, arctgt, 2, GMP_RNDN); /* exact */
+          realprec = prec - 2;
+          goto can_round;
+        }
+
       /* Assignation  */
       MPFR_SET_ZERO (arctgt);
-      twopoweri = 1<<0;
+      twopoweri = 1 << 0;
       MPFR_ASSERTD (n0 >= 4);
       /* FIXME: further reduce the argument so that it is less than
          1/n where n is the output precision. In such a way, the
@@ -366,6 +380,7 @@ mpfr_atan (mpfr_ptr atan, mpfr_srcptr x, mp_rnd_t rnd_mode)
         }
       MPFR_SET_POS (arctgt);
 
+    can_round:
       if (MPFR_LIKELY (MPFR_CAN_ROUND (arctgt, realprec, MPFR_PREC (atan),
                                        rnd_mode)))
         break;