Add special values processing in mpc_sin as in the C99 standard

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

1 files changed, 69 insertions, 14 deletions

diff --git a/src/sin.c b/src/sin.c
index ad1d6ef..047a35e 100644
--- a/src/sin.c
+++ b/src/sin.c
@@ -31,20 +31,61 @@ mpc_sin (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd)
   mp_prec_t prec;
   int ok = 0;
 
-  /* let op = a + i*b, then sin(op) = sin(a)*cosh(b) + i*cos(a)*sinh(b).
+  /* special values */
+  if (!mpfr_number_p (MPC_RE (op)) || !mpfr_number_p (MPC_IM (op)))
+    {
+      if (mpfr_nan_p (MPC_RE (op)) || mpfr_nan_p (MPC_IM (op)))
+        {
+          mpc_set (rop, op, rnd);
 
-     We use the following algorithm (same for the imaginary part),
-     with rounding to nearest for all operations, and working precision w:
+          if (mpfr_nan_p (MPC_IM (op)))
+            /* OP = x +i NaN, then ROP = NaN +i*NaN except for x=0 */
+            {
+              if (!mpfr_zero_p (MPC_RE (op)))
+                mpfr_set_nan (MPC_RE (rop));
+              return;
+            }
+
+          /* OP = NaN +i*y, then ROP = NaN +i*NaN except if x=0 or infinity */
+          if (!mpfr_inf_p (MPC_IM (op)) && !mpfr_zero_p (MPC_IM (op)))
+            mpfr_set_nan (MPC_IM (rop));
+
+          return;
+        }
+
+      if (mpfr_inf_p (MPC_RE (op)))
+        /* OP = infinity +i*y where y is not NaN */
+        {
+          mpfr_set_nan (MPC_RE (rop));
+
+          if (!mpfr_inf_p (MPC_IM (op)) && !mpfr_zero_p (MPC_IM (op)))
+            mpfr_set_nan (MPC_IM (rop));
+          else
+            mpfr_set (MPC_IM (rop), MPC_IM (op), MPC_RND_IM (rnd));
+
+          return;
+        }
+
+      /* OP = x +i*infinity where x is finite */
+      if (mpfr_zero_p (MPC_RE (op)))
+        {
+          mpc_set (rop, op, rnd);
+
+          return;
+        }
+
+      /* ROP = infinity*(sin x -i*cos x) */
+      mpfr_init (x);
+      mpfr_init (y);
+      mpfr_sin_cos (x, y, MPC_RE (op), GMP_RNDZ);
+      mpfr_set_inf (MPC_RE (rop), mpfr_sgn (x));
+      mpfr_set_inf (MPC_IM (rop), mpfr_sgn (y));
+      mpfr_clear (y);
+      mpfr_clear(x);
+
+      return;
+    }
 
-     (1) x = o(sin(a))
-     (2) y = o(cosh(b))
-     (3) r = o(x*y)     
-     then the error on r is at most 4 ulps, since we can write
-     r = sin(a)*cosh(b)*(1+t)^3 with |t| <= 2^(-w),
-     thus for w >= 2, r = sin(a)*cosh(b)*(1+4*t) with |t| <= 2^(-w),
-     thus the relative error is bounded by 4*2^(-w) <= 4*ulp(r).
-  */
-  
   /* special case when the input is real: sin(x) = sin(x) */
   if (mpfr_cmp_ui (MPC_IM(op), 0) == 0)
     {
@@ -61,6 +102,20 @@ mpc_sin (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd)
       return;
     }
 
+  /* let op = a + i*b, then sin(op) = sin(a)*cosh(b) + i*cos(a)*sinh(b).
+
+     We use the following algorithm (same for the imaginary part),
+     with rounding to nearest for all operations, and working precision w:
+
+     (1) x = o(sin(a))
+     (2) y = o(cosh(b))
+     (3) r = o(x*y)     
+     then the error on r is at most 4 ulps, since we can write
+     r = sin(a)*cosh(b)*(1+t)^3 with |t| <= 2^(-w),
+     thus for w >= 2, r = sin(a)*cosh(b)*(1+4*t) with |t| <= 2^(-w),
+     thus the relative error is bounded by 4*2^(-w) <= 4*ulp(r).
+  */
+  
   prec = MPC_MAX_PREC(rop);
 
   mpfr_init2 (x, 2);
@@ -82,8 +137,8 @@ mpc_sin (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd)
                            MPFR_PREC(MPC_RE(rop)));
       if (ok) /* compute imaginary part */
         {
-	  mpfr_sinh (z, MPC_IM(op), GMP_RNDN);
-	  mpfr_mul (y, y, z, GMP_RNDN);
+          mpfr_sinh (z, MPC_IM(op), GMP_RNDN);
+          mpfr_mul (y, y, z, GMP_RNDN);
           ok = mpfr_can_round (y, prec - 2, GMP_RNDN, MPC_RND_IM(rnd),
                                MPFR_PREC(MPC_IM(rop)));
         }