sin_cos: bug fix to allow overlap in argument and result

   previously, there was a problem when computing sin and cos (none of
   rop_sin and rop_cos being NULL), the arguments were non-finite, purely
   real or purely imaginary, and op was equal to one of the output values
tests: added generic and reuse tests for functions with one input and two
   outputs, labelled by "cc_c"; in particular, new test file tsin_cos.c


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

1 files changed, 103 insertions, 72 deletions

diff --git a/src/sin_cos.c b/src/sin_cos.c
index 8b965fb..f09b173 100644
--- a/src/sin_cos.c
+++ b/src/sin_cos.c
@@ -26,16 +26,27 @@ mpc_sin_cos_nonfinite (mpc_ptr rop_sin, mpc_ptr rop_cos, mpc_srcptr op,
    mpc_rnd_t rnd_sin, mpc_rnd_t rnd_cos)
    /* assumes that op (that is, its real or imaginary part) is not finite */
 {
+   int overlap;
+   mpc_t op_loc;
+
+   overlap = (rop_sin == op || rop_cos == op);
+   if (overlap) {
+      mpc_init3 (op_loc, MPC_PREC_RE (op), MPC_PREC_IM (op));
+      mpc_set (op_loc, op, MPC_RNDNN);
+   }
+   else
+      op_loc [0] = op [0];
+
    if (rop_sin != NULL) {
-      if (mpfr_nan_p (MPC_RE (op)) || mpfr_nan_p (MPC_IM (op))) {
-         mpc_set (rop_sin, op, rnd_sin);
-         if (mpfr_nan_p (MPC_IM (op))) {
+      if (mpfr_nan_p (MPC_RE (op_loc)) || mpfr_nan_p (MPC_IM (op_loc))) {
+         mpc_set (rop_sin, op_loc, rnd_sin);
+         if (mpfr_nan_p (MPC_IM (op_loc))) {
             /* sin(x +i*NaN) = NaN +i*NaN, except for x=0 */
             /* sin(-0 +i*NaN) = -0 +i*NaN */
             /* sin(+0 +i*NaN) = +0 +i*NaN */
-            if (!mpfr_zero_p (MPC_RE (op)))
+            if (!mpfr_zero_p (MPC_RE (op_loc)))
                mpfr_set_nan (MPC_RE (rop_sin));
-            else if (!mpfr_inf_p (MPC_IM (op)) && !mpfr_zero_p (MPC_IM (op)))
+            else if (!mpfr_inf_p (MPC_IM (op_loc)) && !mpfr_zero_p (MPC_IM (op_loc)))
                /* sin(NaN -i*Inf) = NaN -i*Inf */
                /* sin(NaN -i*0) = NaN -i*0 */
                /* sin(NaN +i*0) = NaN +i*0 */
@@ -44,10 +55,10 @@ mpc_sin_cos_nonfinite (mpc_ptr rop_sin, mpc_ptr rop_cos, mpc_srcptr op,
                mpfr_set_nan (MPC_IM (rop_sin));
             }
          }
-      else if (mpfr_inf_p (MPC_RE (op))) {
+      else if (mpfr_inf_p (MPC_RE (op_loc))) {
          mpfr_set_nan (MPC_RE (rop_sin));
 
-         if (!mpfr_inf_p (MPC_IM (op)) && !mpfr_zero_p (MPC_IM (op)))
+         if (!mpfr_inf_p (MPC_IM (op_loc)) && !mpfr_zero_p (MPC_IM (op_loc)))
             /* sin(+/-Inf -i*Inf) = NaN -i*Inf */
             /* sin(+/-Inf +i*Inf) = NaN +i*Inf */
             /* sin(+/-Inf +i*y) = NaN +i*NaN, when 0<|y|<Inf */
@@ -55,14 +66,14 @@ mpc_sin_cos_nonfinite (mpc_ptr rop_sin, mpc_ptr rop_cos, mpc_srcptr op,
          else
             /* sin(+/-Inf -i*0) = NaN -i*0 */
             /* sin(+/-Inf +i*0) = NaN +i*0 */
-            mpfr_set (MPC_IM (rop_sin), MPC_IM (op), MPC_RND_IM (rnd_sin));
+            mpfr_set (MPC_IM (rop_sin), MPC_IM (op_loc), MPC_RND_IM (rnd_sin));
       }
-      else if (mpfr_zero_p (MPC_RE (op))) {
+      else if (mpfr_zero_p (MPC_RE (op_loc))) {
          /* sin(-0 -i*Inf) = -0 -i*Inf */
          /* sin(+0 -i*Inf) = +0 -i*Inf */
          /* sin(-0 +i*Inf) = -0 +i*Inf */
          /* sin(+0 +i*Inf) = +0 +i*Inf */
-         mpc_set (rop_sin, op, rnd_sin);
+         mpc_set (rop_sin, op_loc, rnd_sin);
       }
       else {
          /* sin(x -i*Inf) = +Inf*(sin(x) -i*cos(x)) */
@@ -70,46 +81,46 @@ mpc_sin_cos_nonfinite (mpc_ptr rop_sin, mpc_ptr rop_cos, mpc_srcptr op,
          mpfr_t s, c;
          mpfr_init2 (s, 2);
          mpfr_init2 (c, 2);
-         mpfr_sin_cos (s, c, MPC_RE (op), GMP_RNDZ);
+         mpfr_sin_cos (s, c, MPC_RE (op_loc), GMP_RNDZ);
          mpfr_set_inf (MPC_RE (rop_sin), MPFR_SIGN (s));
-         mpfr_set_inf (MPC_IM (rop_sin), MPFR_SIGN (c)*MPFR_SIGN (MPC_IM (op)));
+         mpfr_set_inf (MPC_IM (rop_sin), MPFR_SIGN (c)*MPFR_SIGN (MPC_IM (op_loc)));
          mpfr_clear (s);
          mpfr_clear (c);
       }
    }
 
    if (rop_cos != NULL) {
-      if (mpfr_nan_p (MPC_RE (op))) {
+      if (mpfr_nan_p (MPC_RE (op_loc))) {
          /* cos(NaN + i * NaN) = NaN + i * NaN */
          /* cos(NaN - i * Inf) = +Inf + i * NaN */
          /* cos(NaN + i * Inf) = +Inf + i * NaN */
          /* cos(NaN - i * 0) = NaN - i * 0 */
          /* cos(NaN + i * 0) = NaN + i * 0 */
          /* cos(NaN + i * y) = NaN + i * NaN, when y != 0 */
-         if (mpfr_inf_p (MPC_IM (op)))
+         if (mpfr_inf_p (MPC_IM (op_loc)))
             mpfr_set_inf (MPC_RE (rop_cos), +1);
          else
             mpfr_set_nan (MPC_RE (rop_cos));
 
-         if (mpfr_zero_p (MPC_IM (op)))
-            mpfr_set (MPC_IM (rop_cos), MPC_IM (op), MPC_RND_IM (rnd_cos));
+         if (mpfr_zero_p (MPC_IM (op_loc)))
+            mpfr_set (MPC_IM (rop_cos), MPC_IM (op_loc), MPC_RND_IM (rnd_cos));
          else
             mpfr_set_nan (MPC_IM (rop_cos));
       }
-      else if (mpfr_nan_p (MPC_IM (op))) {
+      else if (mpfr_nan_p (MPC_IM (op_loc))) {
           /* cos(-Inf + i * NaN) = NaN + i * NaN */
           /* cos(+Inf + i * NaN) = NaN + i * NaN */
           /* cos(-0 + i * NaN) = NaN - i * 0 */
           /* cos(+0 + i * NaN) = NaN + i * 0 */
           /* cos(x + i * NaN) = NaN + i * NaN, when x != 0 */
-         if (mpfr_zero_p (MPC_RE (op)))
-            mpfr_set (MPC_IM (rop_cos), MPC_RE (op), MPC_RND_IM (rnd_cos));
+         if (mpfr_zero_p (MPC_RE (op_loc)))
+            mpfr_set (MPC_IM (rop_cos), MPC_RE (op_loc), MPC_RND_IM (rnd_cos));
          else
             mpfr_set_nan (MPC_IM (rop_cos));
 
          mpfr_set_nan (MPC_RE (rop_cos));
       }
-      else if (mpfr_inf_p (MPC_RE (op))) {
+      else if (mpfr_inf_p (MPC_RE (op_loc))) {
          /* cos(-Inf -i*Inf) = cos(+Inf +i*Inf) = -Inf +i*NaN */
          /* cos(-Inf +i*Inf) = cos(+Inf -i*Inf) = +Inf +i*NaN */
          /* cos(-Inf -i*0) = cos(+Inf +i*0) = NaN -i*0 */
@@ -119,26 +130,26 @@ mpc_sin_cos_nonfinite (mpc_ptr rop_sin, mpc_ptr rop_cos, mpc_srcptr op,
          /* SAME_SIGN is useful only if op == (+/-)Inf + i * (+/-)0, but, as
             MPC_RE(OP) might be erased (when ROP == OP), we compute it now */
          const int same_sign =
-            mpfr_signbit (MPC_RE (op)) == mpfr_signbit (MPC_IM (op));
+            mpfr_signbit (MPC_RE (op_loc)) == mpfr_signbit (MPC_IM (op_loc));
 
-         if (mpfr_inf_p (MPC_IM (op)))
+         if (mpfr_inf_p (MPC_IM (op_loc)))
             mpfr_set_inf (MPC_RE (rop_cos), (same_sign ? -1 : +1));
          else
             mpfr_set_nan (MPC_RE (rop_cos));
 
-         if (mpfr_zero_p (MPC_IM (op)))
-            mpfr_setsign (MPC_IM (rop_cos), MPC_IM (op), same_sign,
+         if (mpfr_zero_p (MPC_IM (op_loc)))
+            mpfr_setsign (MPC_IM (rop_cos), MPC_IM (op_loc), same_sign,
                           MPC_RND_IM(rnd_cos));
          else
             mpfr_set_nan (MPC_IM (rop_cos));
       }
-      else if (mpfr_zero_p (MPC_RE (op))) {
+      else if (mpfr_zero_p (MPC_RE (op_loc))) {
          /* cos(-0 -i*Inf) = cos(+0 +i*Inf) = +Inf -i*0 */
          /* cos(-0 +i*Inf) = cos(+0 -i*Inf) = +Inf +i*0 */
          const int same_sign =
-            mpfr_signbit (MPC_RE (op)) == mpfr_signbit (MPC_IM (op));
+            mpfr_signbit (MPC_RE (op_loc)) == mpfr_signbit (MPC_IM (op_loc));
 
-         mpfr_setsign (MPC_IM (rop_cos), MPC_RE (op), same_sign,
+         mpfr_setsign (MPC_IM (rop_cos), MPC_RE (op_loc), same_sign,
                        MPC_RND_IM (rnd_cos));
          mpfr_set_inf (MPC_RE (rop_cos), +1);
       }
@@ -148,15 +159,18 @@ mpc_sin_cos_nonfinite (mpc_ptr rop_sin, mpc_ptr rop_cos, mpc_srcptr op,
          mpfr_t s, c;
          mpfr_init2 (c, 2);
          mpfr_init2 (s, 2);
-         mpfr_sin_cos (s, c, MPC_RE (op), GMP_RNDN);
+         mpfr_sin_cos (s, c, MPC_RE (op_loc), GMP_RNDN);
          mpfr_set_inf (MPC_RE (rop_cos), mpfr_sgn (c));
          mpfr_set_inf (MPC_IM (rop_cos),
-            (mpfr_sgn (MPC_IM (op)) == mpfr_sgn (s) ? -1 : +1));
+            (mpfr_sgn (MPC_IM (op_loc)) == mpfr_sgn (s) ? -1 : +1));
          mpfr_clear (s);
          mpfr_clear (c);
       }
    }
 
+   if (overlap)
+      mpc_clear (op_loc);
+
    return MPC_INEX12 (MPC_INEX (0,0), MPC_INEX (0,0));
       /* everything is exact */
 }
@@ -167,45 +181,56 @@ mpc_sin_cos_real (mpc_ptr rop_sin, mpc_ptr rop_cos, mpc_srcptr op,
    mpc_rnd_t rnd_sin, mpc_rnd_t rnd_cos)
    /* assumes that op is real */
 {
-   int inex_sin_re, inex_cos_re;
+   int inex_sin_re = 0, inex_cos_re = 0;
+      /* Until further notice, assume computations exact; in particular,
+         by definition, for not computed values.                         */
+   mpfr_t s, c;
+   int inex_s, inex_c;
+   int sign_im_op = mpfr_signbit (MPC_IM (op));
+
+   /* sin(x +-0*i) = sin(x) +-0*i*sign(cos(x)) */
+   /* cos(x +-i*0) = cos(x) -+i*0*sign(sin(x)) */
+   if (rop_sin != 0)
+      mpfr_init2 (s, MPC_PREC_RE (rop_sin));
+   else
+      mpfr_init2 (s, 2); /* We need only the sign. */
+   if (rop_cos != NULL)
+      mpfr_init2 (c, MPC_PREC_RE (rop_cos));
+   else
+      mpfr_init2 (c, 2);
+   inex_s = mpfr_sin (s, MPC_RE (op), MPC_RND_RE (rnd_sin));
+   inex_c = mpfr_cos (c, MPC_RE (op), MPC_RND_RE (rnd_cos));
+      /* We cannot use mpfr_sin_cos since we may need two distinct rounding
+         modes and the exact return values. If we need only the sign, an
+         arbitrary rounding mode will work.                                 */
 
    if (rop_sin != NULL) {
-      /* sin(x +-0*i) = sin(x) +-0*i*sign(cos(x)) */
-      mpfr_t (sign);
-      mpfr_init2 (sign, 2);
-      mpfr_cos (sign, MPC_RE (op), GMP_RNDN);
-      if (mpfr_signbit (MPC_IM (op)))
-         MPFR_CHANGE_SIGN (sign);
-
-      inex_sin_re = mpfr_sin (MPC_RE (rop_sin), MPC_RE (op), MPC_RND_RE (rnd_sin));
-
+      mpfr_set (MPC_RE (rop_sin), s, GMP_RNDN); /* exact */
+      inex_sin_re = inex_s;
       mpfr_set_ui (MPC_IM (rop_sin), 0ul, GMP_RNDN);
-      if (mpfr_signbit (sign))
+      if (   ( sign_im_op && !mpfr_signbit (c))
+          || (!sign_im_op &&  mpfr_signbit (c)))
          MPFR_CHANGE_SIGN (MPC_IM (rop_sin));
-
-      mpfr_clear (sign);
+      /* FIXME: simpler implementation with mpfr-3:
+         mpfr_set_zero (MPC_IM (rop_sin),
+            (   ( mpfr_signbit (MPC_IM(op)) && !mpfr_signbit(c))
+             || (!mpfr_signbit (MPC_IM(op)) &&  mpfr_signbit(c)) ? -1 : 1);
+         there is no need to use the variable sign_im_op then, needed now in
+         the case rop_sin == op                                              */
    }
-   else
-      inex_sin_re = 0;
 
    if (rop_cos != NULL) {
-      /* cos(x +-i*0) = cos(x) -+i*0*sign(sin(x)) */
-      mpfr_t sign;
-      mpfr_init2 (sign, 2);
-      mpfr_sin (sign, MPC_RE (op), GMP_RNDN);
-      if (!mpfr_signbit (MPC_IM (op)))
-         MPFR_CHANGE_SIGN (sign);
-
-      inex_cos_re = mpfr_cos (MPC_RE (rop_cos), MPC_RE (op), MPC_RND_RE (rnd_cos));
-
+      mpfr_set (MPC_RE (rop_cos), c, GMP_RNDN); /* exact */
+      inex_cos_re = inex_c;
       mpfr_set_ui (MPC_IM (rop_cos), 0ul, GMP_RNDN);
-      if (mpfr_signbit (sign))
+      if (   ( sign_im_op &&  mpfr_signbit (s))
+          || (!sign_im_op && !mpfr_signbit (s)))
          MPFR_CHANGE_SIGN (MPC_IM (rop_cos));
-
-      mpfr_clear (sign);
+      /* FIXME: see previous MPFR_CHANGE_SIGN */
    }
-   else
-      inex_cos_re = 0;
+
+   mpfr_clear (s);
+   mpfr_clear (c);
 
    return MPC_INEX12 (MPC_INEX (inex_sin_re, 0), MPC_INEX (inex_cos_re, 0));
 }
@@ -216,36 +241,42 @@ mpc_sin_cos_imag (mpc_ptr rop_sin, mpc_ptr rop_cos, mpc_srcptr op,
    mpc_rnd_t rnd_sin, mpc_rnd_t rnd_cos)
    /* assumes that op is purely imaginary */
 {
-   int inex_sin_im, inex_cos_re;
+   int inex_sin_im = 0, inex_cos_re = 0;
+      /* assume exact if not computed */
+   int overlap;
+   mpc_t op_loc;
+
+   overlap = (rop_sin == op || rop_cos == op);
+   if (overlap) {
+      mpc_init3 (op_loc, MPC_PREC_RE (op), MPC_PREC_IM (op));
+      mpc_set (op_loc, op, MPC_RNDNN);
+   }
+   else
+      op_loc [0] = op [0];
 
    if (rop_sin != NULL) {
       /* sin(+-O +i*y) = +-0 +i*sinh(y) */
-      mpfr_set (MPC_RE(rop_sin), MPC_RE(op), GMP_RNDN);
-      inex_sin_im = mpfr_sinh (MPC_IM(rop_sin), MPC_IM(op), MPC_RND_IM(rnd_sin));
+      mpfr_set (MPC_RE(rop_sin), MPC_RE(op_loc), GMP_RNDN);
+      inex_sin_im = mpfr_sinh (MPC_IM(rop_sin), MPC_IM(op_loc), MPC_RND_IM(rnd_sin));
    }
-   else
-      inex_sin_im = 0;
 
    if (rop_cos != NULL) {
       /* cos(-0 - i * y) = cos(+0 + i * y) = cosh(y) - i * 0,
          cos(-0 + i * y) = cos(+0 - i * y) = cosh(y) + i * 0,
          where y >= 0 */
 
-      /* When rop == op, the sign of 0 will be erased, so use it now. */
-      const int imag_sign =
-         mpfr_signbit (MPC_RE (op)) ==  mpfr_signbit (MPC_IM (op));
-
-      if (mpfr_zero_p (MPC_IM (op)))
+      if (mpfr_zero_p (MPC_IM (op_loc)))
         inex_cos_re = mpfr_set_ui (MPC_RE (rop_cos), 1ul, MPC_RND_RE (rnd_cos));
       else
-        inex_cos_re = mpfr_cosh (MPC_RE (rop_cos), MPC_IM (op), MPC_RND_RE (rnd_cos));
+        inex_cos_re = mpfr_cosh (MPC_RE (rop_cos), MPC_IM (op_loc), MPC_RND_RE (rnd_cos));
 
       mpfr_set_ui (MPC_IM (rop_cos), 0ul, MPC_RND_IM (rnd_cos));
-      if (imag_sign)
+      if (mpfr_signbit (MPC_RE (op_loc)) ==  mpfr_signbit (MPC_IM (op_loc)))
          MPFR_CHANGE_SIGN (MPC_IM (rop_cos));
    }
-   else
-      inex_cos_re = 0;
+
+   if (overlap)
+      mpc_clear (op_loc);
 
    return MPC_INEX12 (MPC_INEX (0, inex_sin_im), MPC_INEX (inex_cos_re, 0));
 }