fixed integer undefined behaviors reported by John Regehr (#10838)

git-svn-id: svn://scm.gforge.inria.fr/svn/mpc/trunk@817 211d60ee-9f03-0410-a15a-8952a2c7a4e4

5 files changed, 80 insertions, 42 deletions

diff --git a/BUGS b/BUGS
index f4e260d..8e1277f 100644
--- a/BUGS
+++ b/BUGS
@@ -9,7 +9,3 @@
   sin(x)*sinh(y) is representable. If furthermore an underflow occurred
   in sin(x) (which has not been observed in practice), then the return
   value would be NaN*(+-inf)=NaN.
-  As another example, tan is computed as sin/cos; if there is an overflow
-  in both sin and cos, then NaN+i*NaN is returned even if the result
-  may be representable.
-
diff --git a/src/acos.c b/src/acos.c
index 416eeab..c5f5afe 100644
--- a/src/acos.c
+++ b/src/acos.c
@@ -1,6 +1,6 @@
 /* mpc_acos -- arccosine of a complex number.
 
-Copyright (C) 2009 Philippe Th\'eveny, Paul Zimmermann
+Copyright (C) 2009, 2010 Philippe Th\'eveny, Paul Zimmermann
 
 This file is part of the MPC Library.
 
@@ -206,19 +206,22 @@ mpc_acos (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd)
       inex_im = MPC_INEX_IM(inex); /* inex_im is in {-1, 0, 1} */
       e2 = mpfr_get_exp (MPC_RE(z1));
       mpfr_sub (MPC_RE(z1), pi_over_2, MPC_RE(z1), GMP_RNDN);
-      /* the error on x=Re(z1) is bounded by 1/2 ulp(x) + 2^(e1-p-1) +
-         2^(e2-p-1) */
-      e1 = e1 >= e2 ? e1 + 1 : e2 + 1;
-      /* the error on x is bounded by 1/2 ulp(x) + 2^(e1-p-1) */
-      e1 -= mpfr_get_exp (MPC_RE(z1));
-      /* the error on x is bounded by 1/2 ulp(x) [1 + 2^e1] */
-      e1 = e1 <= 0 ? 0 : e1;
-      /* the error on x is bounded by 2^e1 * ulp(x) */
-      mpfr_neg (MPC_IM(z1), MPC_IM(z1), GMP_RNDN); /* exact */
-      inex_im = -inex_im;
-      if (mpfr_can_round (MPC_RE(z1), p - e1, GMP_RNDN, GMP_RNDZ,
-                          p_re + (MPC_RND_RE(rnd) == GMP_RNDN)))
-        break;
+      if (!mpfr_zero_p (MPC_RE(z1)))
+        {
+          /* the error on x=Re(z1) is bounded by 1/2 ulp(x) + 2^(e1-p-1) +
+             2^(e2-p-1) */
+          e1 = e1 >= e2 ? e1 + 1 : e2 + 1;
+          /* the error on x is bounded by 1/2 ulp(x) + 2^(e1-p-1) */
+          e1 -= mpfr_get_exp (MPC_RE(z1));
+          /* the error on x is bounded by 1/2 ulp(x) [1 + 2^e1] */
+          e1 = e1 <= 0 ? 0 : e1;
+          /* the error on x is bounded by 2^e1 * ulp(x) */
+          mpfr_neg (MPC_IM(z1), MPC_IM(z1), GMP_RNDN); /* exact */
+          inex_im = -inex_im;
+          if (mpfr_can_round (MPC_RE(z1), p - e1, GMP_RNDN, GMP_RNDZ,
+                              p_re + (MPC_RND_RE(rnd) == GMP_RNDN)))
+            break;
+        }
     }
   inex = mpc_set (rop, z1, rnd);
   inex_re = MPC_INEX_RE(inex);
diff --git a/src/atan.c b/src/atan.c
index b393fc2..feb8c47 100644
--- a/src/atan.c
+++ b/src/atan.c
@@ -1,6 +1,6 @@
 /* mpc_atan -- arctangent of a complex number.
 
-Copyright (C) 2009 Philippe Th\'eveny, Paul Zimmermann
+Copyright (C) 2009, 2010 Philippe Th\'eveny, Paul Zimmermann
 
 This file is part of the MPC Library.
 
@@ -331,21 +331,26 @@ mpc_atan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd)
 
         mpfr_sub (y, a, b, GMP_RNDU);
 
-        expo = MPC_MAX (mpfr_get_exp (a), mpfr_get_exp (b));
-        expo = expo - mpfr_get_exp (y) + 1;
-        err = 3 - mpfr_get_exp (y);
-        /* error(j) <= [1 + 2^expo + 7*2^err] ulp(j) */
-        if (expo <= err) /* error(j) <= [1 + 2^{err+1}] ulp(j) */
-          err = (err < 0) ? 1 : err + 2;
-        else
-          err = (expo < 0) ? 1 : expo + 2;
-
-        mpfr_div_2ui (y, y, 2, GMP_RNDN);
         if (mpfr_zero_p (y))
-          err = 2; /* underflow */
-
-        ok = mpfr_can_round (y, p - err, GMP_RNDU, GMP_RNDD,
-                             prec + (MPC_RND_IM (rnd) == GMP_RNDN));
+          ok = 0;
+        else
+          {
+            expo = MPC_MAX (mpfr_get_exp (a), mpfr_get_exp (b));
+            expo = expo - mpfr_get_exp (y) + 1;
+            err = 3 - mpfr_get_exp (y);
+            /* error(j) <= [1 + 2^expo + 7*2^err] ulp(j) */
+            if (expo <= err) /* error(j) <= [1 + 2^{err+1}] ulp(j) */
+              err = (err < 0) ? 1 : err + 2;
+            else
+              err = (expo < 0) ? 1 : expo + 2;
+
+            mpfr_div_2ui (y, y, 2, GMP_RNDN);
+            if (mpfr_zero_p (y))
+              err = 2; /* underflow */
+
+            ok = mpfr_can_round (y, p - err, GMP_RNDU, GMP_RNDD,
+                                 prec + (MPC_RND_IM (rnd) == GMP_RNDN));
+          }
       } while (ok == 0);
 
     inex = mpc_set_fr_fr (rop, x, y, rnd);
diff --git a/src/tan.c b/src/tan.c
index 1f2714c..da8b70c 100644
--- a/src/tan.c
+++ b/src/tan.c
@@ -1,6 +1,6 @@
 /* mpc_tan -- tangent of a complex number.
 
-Copyright (C) 2008, 2009, 2010 Philippe Th\'eveny, Andreas Enge
+Copyright (C) 2008, 2009, 2010 Philippe Th\'eveny, Andreas Enge, Paul Zimmermann
 
 This file is part of the MPC Library.
 
@@ -19,6 +19,8 @@ along with the MPC 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. */
 
+#include <stdio.h>    /* for MPC_ASSERT */
+#include <limits.h>
 #include "mpc-impl.h"
 
 int
@@ -181,6 +183,7 @@ mpc_tan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd)
   do
     {
       mpfr_exp_t k, exr, eyr, eyi, ezr;
+      int isinf = 0;
 
       ok = 0;
 
@@ -193,16 +196,49 @@ mpc_tan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd)
          above), sin x and cos y are never exact, so rounding away from 0 is
          rounding towards 0 and adding one ulp to the absolute value */
       mpc_sin (x, op, MPC_RNDZZ);
+      isinf = mpfr_inf_p (MPC_RE (x)) || mpfr_inf_p (MPC_IM (x));
+      /* if the real or imaginary part of x is infinite, it means that Im(op)
+         was large, in which case the result is
+         sign(tan(Re(op)))*0 + sign(Im(op))*I,
+         where sign(tan(Re(op))) = sign(Re(x))*sign(Re(y)) */
       mpfr_signbit (MPC_RE (x)) ?
         mpfr_nextbelow (MPC_RE (x)) : mpfr_nextabove (MPC_RE (x));
       mpfr_signbit (MPC_IM (x)) ?
         mpfr_nextbelow (MPC_IM (x)) : mpfr_nextabove (MPC_IM (x));
+      MPC_ASSERT (mpfr_zero_p (MPC_RE (x)) == 0);
       exr = mpfr_get_exp (MPC_RE (x));
       mpc_cos (y, op, MPC_RNDZZ);
+      isinf = isinf || mpfr_inf_p (MPC_RE (y)) || mpfr_inf_p (MPC_IM (y));
       mpfr_signbit (MPC_RE (y)) ?
         mpfr_nextbelow (MPC_RE (y)) : mpfr_nextabove (MPC_RE (y));
       mpfr_signbit (MPC_IM (y)) ?
         mpfr_nextbelow (MPC_IM (y)) : mpfr_nextabove (MPC_IM (y));
+
+      if (isinf)
+        {
+          int inex_re, inex_im;
+          mpfr_set_ui (MPC_RE (rop), 0, GMP_RNDN);
+          if (mpfr_sgn (MPC_RE (x)) * mpfr_sgn (MPC_RE (y)) < 0)
+            {
+              mpfr_neg (MPC_RE (rop), MPC_RE (rop), GMP_RNDN);
+              inex_re = 1;
+            }
+          else
+            inex_re = -1; /* +0 is rounded down */
+          if (mpfr_sgn (MPC_IM (op)) > 0)
+            {
+              mpfr_set_ui (MPC_IM (rop), 1, GMP_RNDN);
+              inex_im = 1;
+            }
+          else
+            {
+              mpfr_set_si (MPC_IM (rop), -1, GMP_RNDN);
+              inex_im = -1;
+            }
+          inex = MPC_INEX(inex_re, inex_im);
+          goto end;
+        }
+
       eyr = mpfr_get_exp (MPC_RE (y));
       eyi = mpfr_get_exp (MPC_IM (y));
 
@@ -225,6 +261,7 @@ mpc_tan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd)
       if (MPC_INEX_IM (inex))
         mpfr_signbit (MPC_IM (x)) ?
           mpfr_nextbelow (MPC_IM (x)) : mpfr_nextabove (MPC_IM (x));
+      MPC_ASSERT (mpfr_zero_p (MPC_RE (x)) == 0);
       ezr = mpfr_get_exp (MPC_RE (x));
 
       /* FIXME: compute
@@ -233,14 +270,14 @@ mpc_tan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd)
       k = exr - ezr + MPC_MAX(-eyr, eyr - 2 * eyi);
       err = k < 2 ? 7 : (k == 2 ? 8 : (5 + k));
 
-      /* Can the real part be rounded ? */
+      /* Can the real part be rounded? */
       ok = (!mpfr_number_p (MPC_RE (x)))
            || mpfr_can_round (MPC_RE(x), prec - err, GMP_RNDN, GMP_RNDZ,
                       MPC_PREC_RE(rop) + (MPC_RND_RE(rnd) == GMP_RNDN));
 
       if (ok)
         {
-          /* Can the imaginary part be rounded ? */
+          /* Can the imaginary part be rounded? */
           ok = (!mpfr_number_p (MPC_IM (x)))
                || mpfr_can_round (MPC_IM(x), prec - 6, GMP_RNDN, GMP_RNDZ,
                       MPC_PREC_IM(rop) + (MPC_RND_IM(rnd) == GMP_RNDN));
@@ -250,6 +287,7 @@ mpc_tan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd)
 
   inex = mpc_set (rop, x, rnd);
 
+ end:
   mpc_clear (x);
   mpc_clear (y);
 
diff --git a/tests/tan.dat b/tests/tan.dat
index 974c882..0a97e48 100644
--- a/tests/tan.dat
+++ b/tests/tan.dat
@@ -126,9 +126,5 @@
 - + 9 -0x9bp-51 9 -1    9 -0x16dp-8 9 -0x77p-3   N N
 
 # huge values
-# should yield finite result as in commented lines, in particular simple imaginary part;
-# in current implementation, yields NaN due to intermediate overflow in sin and cos
-#- + 53 -0xE72F36A66DDF7p-2765858364927957 53 -1 53 0x4580CBF242683p-3 53 -0x1B3E8A3660D279p-3 N N
-#- - 53 -0xD10502370373Fp-441000018035230 53 1 53 -0x1B3E8A3660D279p-3 53 0x4580CBF242683p-3 N N
-0 0 53 nan 53 nan  53 0x4580CBF242683p-3 53 -0x1B3E8A3660D279p-3 N N
-0 0 53 nan 53 nan  53 -0x1B3E8A3660D279p-3 53 0x4580CBF242683p-3 N N
++ - 53 -0 53 -1  53 0x4580CBF242683p-3 53 -0x1B3E8A3660D279p-3 N N
++ + 53 -0 53 +1  53 -0x1B3E8A3660D279p-3 53 0x4580CBF242683p-3 N N