src/tan.c

   1 /* mpc_tan -- tangent of a complex number.
   2
   3 Copyright (C) 2008, 2009 Philippe Th\'eveny, Andreas Enge
   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_tan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd)
  26 {
  27   mpc_t x, y;
  28   mp_prec_t prec;
  29   mp_exp_t err;
  30   int ok = 0;
  31   int inex;
  32
  33   /* special values */
  34   if (!mpfr_number_p (MPC_RE (op)) || !mpfr_number_p (MPC_IM (op)))
  35     {
  36       if (mpfr_nan_p (MPC_RE (op)))
  37         {
  38           if (mpfr_inf_p (MPC_IM (op)))
  39             /* tan(NaN -i*Inf) = +/-0 -i */
  40             /* tan(NaN +i*Inf) = +/-0 +i */
  41             {
  42               /* exact unless 1 is not in exponent range */
  43               inex = mpc_set_si_si (rop, 0,
  44                                     (MPFR_SIGN (MPC_IM (op)) < 0) ? -1 : +1,
  45                                     rnd);
  46             }
  47           else
  48             /* tan(NaN +i*y) = NaN +i*NaN, when y is finite */
  49             /* tan(NaN +i*NaN) = NaN +i*NaN */
  50             {
  51               mpfr_set_nan (MPC_RE (rop));
  52               mpfr_set_nan (MPC_IM (rop));
  53               inex = MPC_INEX (0, 0); /* always exact */
  54             }
  55         }
  56       else if (mpfr_nan_p (MPC_IM (op)))
  57         {
  58           if (mpfr_cmp_ui (MPC_RE (op), 0) == 0)
  59             /* tan(-0 +i*NaN) = -0 +i*NaN */
  60             /* tan(+0 +i*NaN) = +0 +i*NaN */
  61             {
  62               mpc_set (rop, op, rnd);
  63               inex = MPC_INEX (0, 0); /* always exact */
  64             }
  65           else
  66             /* tan(x +i*NaN) = NaN +i*NaN, when x != 0 */
  67             {
  68               mpfr_set_nan (MPC_RE (rop));
  69               mpfr_set_nan (MPC_IM (rop));
  70               inex = MPC_INEX (0, 0); /* always exact */
  71             }
  72         }
  73       else if (mpfr_inf_p (MPC_RE (op)))
  74         {
  75           if (mpfr_inf_p (MPC_IM (op)))
  76             /* tan(-Inf -i*Inf) = -/+0 -i */
  77             /* tan(-Inf +i*Inf) = -/+0 +i */
  78             /* tan(+Inf -i*Inf) = +/-0 -i */
  79             /* tan(+Inf +i*Inf) = +/-0 +i */
  80             {
  81               const int sign_re = mpfr_signbit (MPC_RE (op));
  82               int inex_im;
  83
  84               mpfr_set_ui (MPC_RE (rop), 0, MPC_RND_RE (rnd));
  85               mpfr_setsign (MPC_RE (rop), MPC_RE (rop), sign_re, GMP_RNDN);
  86
  87               /* exact, unless 1 is not in exponent range */
  88               inex_im = mpfr_set_si (MPC_IM (rop),
  89                                      mpfr_signbit (MPC_IM (op)) ? -1 : +1,
  90                                      MPC_RND_IM (rnd));
  91               inex = MPC_INEX (0, inex_im);
  92             }
  93           else
  94             /* tan(-Inf +i*y) = tan(+Inf +i*y) = NaN +i*NaN, when y is
  95                finite */
  96             {
  97               mpfr_set_nan (MPC_RE (rop));
  98               mpfr_set_nan (MPC_IM (rop));
  99               inex = MPC_INEX (0, 0); /* always exact */
 100             }
 101         }
 102       else
 103         /* tan(x -i*Inf) = +0*sin(x)*cos(x) -i, when x is finite */
 104         /* tan(x +i*Inf) = +0*sin(x)*cos(x) +i, when x is finite */
 105         {
 106           mpfr_t c;
 107           mpfr_t s;
 108           int inex_im;
 109
 110           mpfr_init (c);
 111           mpfr_init (s);
 112
 113           mpfr_sin_cos (s, c, MPC_RE (op), GMP_RNDN);
 114           mpfr_set_ui (MPC_RE (rop), 0, MPC_RND_RE (rnd));
 115           mpfr_setsign (MPC_RE (rop), MPC_RE (rop),
 116                         mpfr_signbit (c) != mpfr_signbit (s), GMP_RNDN);
 117           /* exact, unless 1 is not in exponent range */
 118           inex_im = mpfr_set_si (MPC_IM (rop),
 119                                  (mpfr_signbit (MPC_IM (op)) ? -1 : +1),
 120                                  MPC_RND_IM (rnd));
 121           inex = MPC_INEX (0, inex_im);
 122
 123           mpfr_clear (s);
 124           mpfr_clear (c);
 125         }
 126
 127       return inex;
 128     }
 129
 130   if (mpfr_zero_p (MPC_RE (op)))
 131     /* tan(-0 -i*y) = -0 +i*tanh(y), when y is finite. */
 132     /* tan(+0 +i*y) = +0 +i*tanh(y), when y is finite. */
 133     {
 134       int inex_im;
 135
 136       mpfr_set (MPC_RE (rop), MPC_RE (op), MPC_RND_RE (rnd));
 137       inex_im = mpfr_tanh (MPC_IM (rop), MPC_IM (op), MPC_RND_IM (rnd));
 138
 139       return MPC_INEX (0, inex_im);
 140     }
 141
 142   if (mpfr_zero_p (MPC_IM (op)))
 143     /* tan(x -i*0) = tan(x) -i*0, when x is finite. */
 144     /* tan(x +i*0) = tan(x) +i*0, when x is finite. */
 145     {
 146       int inex_re;
 147
 148       inex_re = mpfr_tan (MPC_RE (rop), MPC_RE (op), MPC_RND_RE (rnd));
 149       mpfr_set (MPC_IM (rop), MPC_IM (op), MPC_RND_IM (rnd));
 150
 151       return MPC_INEX (inex_re, 0);
 152     }
 153
 154   /* ordinary (non-zero) numbers */
 155
 156   /* tan(op) = sin(op) / cos(op).
 157
 158      We use the following algorithm with rounding away from 0 for all
 159      operations, and working precision w:
 160
 161      (1) x = A(sin(op))
 162      (2) y = A(cos(op))
 163      (3) z = A(x/y)
 164
 165      the error on Im(z) is at most 81 ulp,
 166      the error on Re(z) is at most
 167      7 ulp if k < 2,
 168      8 ulp if k = 2,
 169      else 5+k ulp, where
 170      k = Exp(Re(x))+Exp(Re(y))-2min{Exp(Re(y)), Exp(Im(y))}-Exp(Re(x/y))
 171      see proof in algorithms.tex.
 172   */
 173
 174   prec = MPC_MAX_PREC(rop);
 175
 176   mpc_init2 (x, 2);
 177   mpc_init2 (y, 2);
 178
 179   err = 7;
 180
 181   do
 182     {
 183       mp_exp_t k;
 184       mp_exp_t exr, eyr, eyi, ezr;
 185
 186       ok = 0;
 187
 188       /* FIXME: prevent addition overflow */
 189       prec += mpc_ceil_log2 (prec) + err;
 190       mpc_set_prec (x, prec);
 191       mpc_set_prec (y, prec);
 192
 193       /* rounding away from zero: except in the cases x=0 or y=0 (processed
 194          above), sin x and cos y are never exact, so rounding away from 0 is
 195          rounding towards 0 and adding one ulp to the absolute value */
 196       mpc_sin (x, op, MPC_RNDZZ);
 197       mpfr_signbit (MPC_RE (x)) ?
 198         mpfr_nextbelow (MPC_RE (x)) : mpfr_nextabove (MPC_RE (x));
 199       mpfr_signbit (MPC_IM (x)) ?
 200         mpfr_nextbelow (MPC_IM (x)) : mpfr_nextabove (MPC_IM (x));
 201       exr = mpfr_get_exp (MPC_RE (x));
 202       mpc_cos (y, op, MPC_RNDZZ);
 203       mpfr_signbit (MPC_RE (y)) ?
 204         mpfr_nextbelow (MPC_RE (y)) : mpfr_nextabove (MPC_RE (y));
 205       mpfr_signbit (MPC_IM (y)) ?
 206         mpfr_nextbelow (MPC_IM (y)) : mpfr_nextabove (MPC_IM (y));
 207       eyr = mpfr_get_exp (MPC_RE (y));
 208       eyi = mpfr_get_exp (MPC_IM (y));
 209
 210       /* some parts of the quotient may be exact */
 211       inex = mpc_div (x, x, y, MPC_RNDZZ);
 212       /* OP is no pure real nor pure imaginary, so in theory the real and
 213          imaginary parts of its tangent cannot be null. However due to
 214          rouding errors this might happen. Consider for example
 215          tan(1+14*I) = 1.26e-10 + 1.00*I. For small precision sin(op) and
 216          cos(op) differ only by a factor I, thus after mpc_div x = I and
 217          its real part is zero. */
 218       if (mpfr_zero_p (MPC_RE (x)) || mpfr_zero_p (MPC_IM (x)))
 219         {
 220           err = prec; /* double precision */
 221           continue;
 222         }
 223       if (MPC_INEX_RE (inex))
 224         mpfr_signbit (MPC_RE (x)) ?
 225           mpfr_nextbelow (MPC_RE (x)) : mpfr_nextabove (MPC_RE (x));
 226       if (MPC_INEX_IM (inex))
 227         mpfr_signbit (MPC_IM (x)) ?
 228           mpfr_nextbelow (MPC_IM (x)) : mpfr_nextabove (MPC_IM (x));
 229       ezr = mpfr_get_exp (MPC_RE (x));
 230
 231       /* FIXME: compute
 232          k = Exp(Re(x))+Exp(Re(y))-2min{Exp(Re(y)), Exp(Im(y))}-Exp(Re(x/y))
 233          avoiding overflow */
 234       k = exr - ezr + MPC_MAX(-eyr, eyr - 2 * eyi);
 235       err = k < 2 ? 7 : (k == 2 ? 8 : (5 + k));
 236
 237       /* Can the real part be rounded ? */
 238       ok = mpfr_inf_p (MPC_RE (x))
 239         || mpfr_can_round (MPC_RE(x), prec - err, GMP_RNDN, GMP_RNDZ,
 240                       MPFR_PREC(MPC_RE(rop)) + (MPC_RND_RE(rnd) == GMP_RNDN));
 241
 242       if (ok)
 243         {
 244           /* Can the imaginary part be rounded ? */
 245           ok = mpfr_inf_p (MPC_IM (x))
 246             || mpfr_can_round (MPC_IM(x), prec - 6, GMP_RNDN, GMP_RNDZ,
 247                       MPFR_PREC(MPC_IM(rop)) + (MPC_RND_IM(rnd) == GMP_RNDN));
 248         }
 249     }
 250   while (ok == 0);
 251
 252   inex = mpc_set (rop, x, rnd);
 253
 254   mpc_clear (x);
 255   mpc_clear (y);
 256
 257   return inex;
 258 }