src/sin.c

   1 /* mpc_sin -- sine of a complex number.
   2
   3 Copyright (C) 2007, 2009 Paul Zimmermann, Philippe Th\'eveny
   4
   5 This file is part of the MPC Library.
   6
   7 The MPC Library is free software; you can redistribute it and/or modify
   8 it under the terms of the GNU Lesser General Public License as published by
   9 the Free Software Foundation; either version 2.1 of the License, or (at your
  10 option) any later version.
  11
  12 The MPC Library is distributed in the hope that it will be useful, but
  13 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  14 or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General Public
  15 License for more details.
  16
  17 You should have received a copy of the GNU Lesser General Public License
  18 along with the MPC Library; see the file COPYING.LIB.  If not, write to
  19 the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
  20 MA 02111-1307, USA. */
  21
  22 #include "mpc-impl.h"
  23
  24 int
  25 mpc_sin (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd)
  26 {
  27   mpfr_t x, y, z;
  28   mp_prec_t prec;
  29   int ok = 0;
  30   int inex_re, inex_im;
  31
  32   /* special values */
  33   if (!mpfr_number_p (MPC_RE (op)) || !mpfr_number_p (MPC_IM (op)))
  34     {
  35       if (mpfr_nan_p (MPC_RE (op)) || mpfr_nan_p (MPC_IM (op)))
  36         {
  37           mpc_set (rop, op, rnd);
  38
  39           if (mpfr_nan_p (MPC_IM (op)))
  40             {
  41               /* sin(x +i*NaN) = NaN +i*NaN, except for x=0 */
  42               /* sin(-0 +i*NaN) = -0 +i*NaN */
  43               /* sin(+0 +i*NaN) = +0 +i*NaN */
  44               if (!mpfr_zero_p (MPC_RE (op)))
  45                 mpfr_set_nan (MPC_RE (rop));
  46               else if (!mpfr_inf_p (MPC_IM (op))
  47                        && !mpfr_zero_p (MPC_IM (op)))
  48                 /* sin(NaN -i*Inf) = NaN -i*Inf */
  49                 /* sin(NaN -i*0) = NaN -i*0 */
  50                 /* sin(NaN +i*0) = NaN +i*0 */
  51                 /* sin(NaN +i*Inf) = NaN +i*Inf */
  52                 /* sin(NaN +i*y) = NaN +i*NaN, when 0<|y|<Inf */
  53                 mpfr_set_nan (MPC_IM (rop));
  54             }
  55         }
  56       else if (mpfr_inf_p (MPC_RE (op)))
  57         {
  58           mpfr_set_nan (MPC_RE (rop));
  59
  60           if (!mpfr_inf_p (MPC_IM (op)) && !mpfr_zero_p (MPC_IM (op)))
  61             /* sin(+/-Inf -i*Inf) = NaN -i*Inf */
  62             /* sin(+/-Inf +i*Inf) = NaN +i*Inf */
  63             /* sin(+/-Inf +i*y) = NaN +i*NaN, when 0<|y|<Inf */
  64             mpfr_set_nan (MPC_IM (rop));
  65           else
  66             /* sin(+/-Inf -i*0) = NaN -i*0 */
  67             /* sin(+/-Inf +i*0) = NaN +i*0 */
  68             mpfr_set (MPC_IM (rop), MPC_IM (op), MPC_RND_IM (rnd));
  69         }
  70       else if (mpfr_zero_p (MPC_RE (op)))
  71         /* sin(-0 -i*Inf) = -0 -i*Inf */
  72         /* sin(+0 -i*Inf) = +0 -i*Inf */
  73         /* sin(-0 +i*Inf) = -0 +i*Inf */
  74         /* sin(+0 +i*Inf) = +0 +i*Inf */
  75         {
  76           mpc_set (rop, op, rnd);
  77         }
  78       else
  79         /* sin(x -i*Inf) = +Inf*(sin(x) -i*cos(x)) */
  80         /* sin(x +i*Inf) = +Inf*(sin(x) +i*cos(x)) */
  81         {
  82           mpfr_init2 (x, 2);
  83           mpfr_init2 (y, 2);
  84           mpfr_sin_cos (x, y, MPC_RE (op), GMP_RNDZ);
  85           mpfr_set_inf (MPC_RE (rop), MPFR_SIGN (x));
  86           mpfr_set_inf (MPC_IM (rop), MPFR_SIGN (y)*MPFR_SIGN (MPC_IM (op)));
  87           mpfr_clear (y);
  88           mpfr_clear(x);
  89         }
  90
  91       return MPC_INEX (0, 0); /* exact in all cases*/
  92     }
  93
  94   /* special case when the input is real: */
  95   /* sin(x -0*i) = sin(x) -0*i*cos(x) */
  96   /* sin(x +0*i) = sin(x) +0*i*cos(x) */
  97   if (mpfr_cmp_ui (MPC_IM(op), 0) == 0)
  98     {
  99       mpfr_init2 (x, 2);
 100       mpfr_cos (x, MPC_RE (op), MPC_RND_RE (rnd));
 101       inex_re = mpfr_sin (MPC_RE (rop), MPC_RE (op), MPC_RND_RE (rnd));
 102       mpfr_mul (MPC_IM(rop), MPC_IM(op), x, MPC_RND_IM(rnd));
 103       mpfr_clear (x);
 104
 105       return MPC_INEX (inex_re, 0);
 106     }
 107
 108   /* special case when the input is imaginary:
 109      sin(+/-O +i*y) = +/-0 +i*sinh(y) */
 110   if (mpfr_cmp_ui (MPC_RE(op), 0) == 0)
 111     {
 112       mpfr_set (MPC_RE(rop), MPC_RE(op), MPC_RND_RE(rnd));
 113       inex_im = mpfr_sinh (MPC_IM(rop), MPC_IM(op), MPC_RND_IM(rnd));
 114
 115       return MPC_INEX (0, inex_im);
 116     }
 117
 118   /* let op = a + i*b, then sin(op) = sin(a)*cosh(b) + i*cos(a)*sinh(b).
 119
 120      We use the following algorithm (same for the imaginary part),
 121      with rounding to nearest for all operations, and working precision w:
 122
 123      (1) x = o(sin(a))
 124      (2) y = o(cosh(b))
 125      (3) r = o(x*y)
 126      then the error on r is at most 4 ulps, since we can write
 127      r = sin(a)*cosh(b)*(1+t)^3 with |t| <= 2^(-w),
 128      thus for w >= 2, r = sin(a)*cosh(b)*(1+4*t) with |t| <= 2^(-w),
 129      thus the relative error is bounded by 4*2^(-w) <= 4*ulp(r).
 130   */
 131
 132   prec = MPC_MAX_PREC(rop);
 133
 134   mpfr_init2 (x, 2);
 135   mpfr_init2 (y, 2);
 136   mpfr_init2 (z, 2);
 137
 138   do
 139     {
 140       prec += mpc_ceil_log2 (prec) + 5;
 141
 142       mpfr_set_prec (x, prec);
 143       mpfr_set_prec (y, prec);
 144       mpfr_set_prec (z, prec);
 145
 146       mpfr_sin_cos (x, y, MPC_RE(op), GMP_RNDN);
 147       mpfr_cosh (z, MPC_IM(op), GMP_RNDN);
 148       mpfr_mul (x, x, z, GMP_RNDN);
 149       ok = mpfr_can_round (x, prec - 2, GMP_RNDN, GMP_RNDZ,
 150                       MPFR_PREC(MPC_RE(rop)) + (MPC_RND_RE(rnd) == GMP_RNDN));
 151       if (ok) /* compute imaginary part */
 152         {
 153           mpfr_sinh (z, MPC_IM(op), GMP_RNDN);
 154           mpfr_mul (y, y, z, GMP_RNDN);
 155           ok = mpfr_can_round (y, prec - 2, GMP_RNDN, GMP_RNDZ,
 156                       MPFR_PREC(MPC_IM(rop)) + (MPC_RND_IM(rnd) == GMP_RNDN));
 157         }
 158     }
 159   while (ok == 0);
 160
 161   inex_re = mpfr_set (MPC_RE(rop), x, MPC_RND_RE(rnd));
 162   inex_im = mpfr_set (MPC_IM(rop), y, MPC_RND_IM(rnd));
 163
 164   mpfr_clear (x);
 165   mpfr_clear (y);
 166   mpfr_clear (z);
 167
 168   return MPC_INEX (inex_re, inex_im);
 169 }