contrib/mpc/src/tan.c

   1 /* mpc_tan -- tangent of a complex number.
   2
   3 Copyright (C) 2008, 2009, 2010, 2011 INRIA
   4
   5 This file is part of GNU MPC.
   6
   7 GNU MPC is free software; you can redistribute it and/or modify it under
   8 the terms of the GNU Lesser General Public License as published by the
   9 Free Software Foundation; either version 3 of the License, or (at your
  10 option) any later version.
  11
  12 GNU MPC is distributed in the hope that it will be useful, but WITHOUT ANY
  13 WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
  14 FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
  15 more details.
  16
  17 You should have received a copy of the GNU Lesser General Public License
  18 along with this program. If not, see http://www.gnu.org/licenses/ .
  19 */
  20
  21 #include <stdio.h>    /* for MPC_ASSERT */
  22 #include <limits.h>
  23 #include "mpc-impl.h"
  24
  25 int
  26 mpc_tan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd)
  27 {
  28   mpc_t x, y;
  29   mpfr_prec_t prec;
  30   mpfr_exp_t err;
  31   int ok = 0;
  32   int inex;
  33
  34   /* special values */
  35   if (!mpc_fin_p (op))
  36     {
  37       if (mpfr_nan_p (mpc_realref (op)))
  38         {
  39           if (mpfr_inf_p (mpc_imagref (op)))
  40             /* tan(NaN -i*Inf) = +/-0 -i */
  41             /* tan(NaN +i*Inf) = +/-0 +i */
  42             {
  43               /* exact unless 1 is not in exponent range */
  44               inex = mpc_set_si_si (rop, 0,
  45                                     (MPFR_SIGN (mpc_imagref (op)) < 0) ? -1 : +1,
  46                                     rnd);
  47             }
  48           else
  49             /* tan(NaN +i*y) = NaN +i*NaN, when y is finite */
  50             /* tan(NaN +i*NaN) = NaN +i*NaN */
  51             {
  52               mpfr_set_nan (mpc_realref (rop));
  53               mpfr_set_nan (mpc_imagref (rop));
  54               inex = MPC_INEX (0, 0); /* always exact */
  55             }
  56         }
  57       else if (mpfr_nan_p (mpc_imagref (op)))
  58         {
  59           if (mpfr_cmp_ui (mpc_realref (op), 0) == 0)
  60             /* tan(-0 +i*NaN) = -0 +i*NaN */
  61             /* tan(+0 +i*NaN) = +0 +i*NaN */
  62             {
  63               mpc_set (rop, op, rnd);
  64               inex = MPC_INEX (0, 0); /* always exact */
  65             }
  66           else
  67             /* tan(x +i*NaN) = NaN +i*NaN, when x != 0 */
  68             {
  69               mpfr_set_nan (mpc_realref (rop));
  70               mpfr_set_nan (mpc_imagref (rop));
  71               inex = MPC_INEX (0, 0); /* always exact */
  72             }
  73         }
  74       else if (mpfr_inf_p (mpc_realref (op)))
  75         {
  76           if (mpfr_inf_p (mpc_imagref (op)))
  77             /* tan(-Inf -i*Inf) = -/+0 -i */
  78             /* tan(-Inf +i*Inf) = -/+0 +i */
  79             /* tan(+Inf -i*Inf) = +/-0 -i */
  80             /* tan(+Inf +i*Inf) = +/-0 +i */
  81             {
  82               const int sign_re = mpfr_signbit (mpc_realref (op));
  83               int inex_im;
  84
  85               mpfr_set_ui (mpc_realref (rop), 0, MPC_RND_RE (rnd));
  86               mpfr_setsign (mpc_realref (rop), mpc_realref (rop), sign_re, GMP_RNDN);
  87
  88               /* exact, unless 1 is not in exponent range */
  89               inex_im = mpfr_set_si (mpc_imagref (rop),
  90                                      mpfr_signbit (mpc_imagref (op)) ? -1 : +1,
  91                                      MPC_RND_IM (rnd));
  92               inex = MPC_INEX (0, inex_im);
  93             }
  94           else
  95             /* tan(-Inf +i*y) = tan(+Inf +i*y) = NaN +i*NaN, when y is
  96                finite */
  97             {
  98               mpfr_set_nan (mpc_realref (rop));
  99               mpfr_set_nan (mpc_imagref (rop));
 100               inex = MPC_INEX (0, 0); /* always exact */
 101             }
 102         }
 103       else
 104         /* tan(x -i*Inf) = +0*sin(x)*cos(x) -i, when x is finite */
 105         /* tan(x +i*Inf) = +0*sin(x)*cos(x) +i, when x is finite */
 106         {
 107           mpfr_t c;
 108           mpfr_t s;
 109           int inex_im;
 110
 111           mpfr_init (c);
 112           mpfr_init (s);
 113
 114           mpfr_sin_cos (s, c, mpc_realref (op), GMP_RNDN);
 115           mpfr_set_ui (mpc_realref (rop), 0, MPC_RND_RE (rnd));
 116           mpfr_setsign (mpc_realref (rop), mpc_realref (rop),
 117                         mpfr_signbit (c) != mpfr_signbit (s), GMP_RNDN);
 118           /* exact, unless 1 is not in exponent range */
 119           inex_im = mpfr_set_si (mpc_imagref (rop),
 120                                  (mpfr_signbit (mpc_imagref (op)) ? -1 : +1),
 121                                  MPC_RND_IM (rnd));
 122           inex = MPC_INEX (0, inex_im);
 123
 124           mpfr_clear (s);
 125           mpfr_clear (c);
 126         }
 127
 128       return inex;
 129     }
 130
 131   if (mpfr_zero_p (mpc_realref (op)))
 132     /* tan(-0 -i*y) = -0 +i*tanh(y), when y is finite. */
 133     /* tan(+0 +i*y) = +0 +i*tanh(y), when y is finite. */
 134     {
 135       int inex_im;
 136
 137       mpfr_set (mpc_realref (rop), mpc_realref (op), MPC_RND_RE (rnd));
 138       inex_im = mpfr_tanh (mpc_imagref (rop), mpc_imagref (op), MPC_RND_IM (rnd));
 139
 140       return MPC_INEX (0, inex_im);
 141     }
 142
 143   if (mpfr_zero_p (mpc_imagref (op)))
 144     /* tan(x -i*0) = tan(x) -i*0, when x is finite. */
 145     /* tan(x +i*0) = tan(x) +i*0, when x is finite. */
 146     {
 147       int inex_re;
 148
 149       inex_re = mpfr_tan (mpc_realref (rop), mpc_realref (op), MPC_RND_RE (rnd));
 150       mpfr_set (mpc_imagref (rop), mpc_imagref (op), MPC_RND_IM (rnd));
 151
 152       return MPC_INEX (inex_re, 0);
 153     }
 154
 155   /* ordinary (non-zero) numbers */
 156
 157   /* tan(op) = sin(op) / cos(op).
 158
 159      We use the following algorithm with rounding away from 0 for all
 160      operations, and working precision w:
 161
 162      (1) x = A(sin(op))
 163      (2) y = A(cos(op))
 164      (3) z = A(x/y)
 165
 166      the error on Im(z) is at most 81 ulp,
 167      the error on Re(z) is at most
 168      7 ulp if k < 2,
 169      8 ulp if k = 2,
 170      else 5+k ulp, where
 171      k = Exp(Re(x))+Exp(Re(y))-2min{Exp(Re(y)), Exp(Im(y))}-Exp(Re(x/y))
 172      see proof in algorithms.tex.
 173   */
 174
 175   prec = MPC_MAX_PREC(rop);
 176
 177   mpc_init2 (x, 2);
 178   mpc_init2 (y, 2);
 179
 180   err = 7;
 181
 182   do
 183     {
 184       mpfr_exp_t k, 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_cos (x, y, op, MPC_RNDZZ, MPC_RNDZZ);
 197       MPFR_ADD_ONE_ULP (mpc_realref (x));
 198       MPFR_ADD_ONE_ULP (mpc_imagref (x));
 199       MPFR_ADD_ONE_ULP (mpc_realref (y));
 200       MPFR_ADD_ONE_ULP (mpc_imagref (y));
 201       MPC_ASSERT (mpfr_zero_p (mpc_realref (x)) == 0);
 202
 203       if (   mpfr_inf_p (mpc_realref (x)) || mpfr_inf_p (mpc_imagref (x))
 204           || mpfr_inf_p (mpc_realref (y)) || mpfr_inf_p (mpc_imagref (y))) {
 205          /* If the real or imaginary part of x is infinite, it means that
 206             Im(op) was large, in which case the result is
 207             sign(tan(Re(op)))*0 + sign(Im(op))*I,
 208             where sign(tan(Re(op))) = sign(Re(x))*sign(Re(y)). */
 209           int inex_re, inex_im;
 210           mpfr_set_ui (mpc_realref (rop), 0, GMP_RNDN);
 211           if (mpfr_sgn (mpc_realref (x)) * mpfr_sgn (mpc_realref (y)) < 0)
 212             {
 213               mpfr_neg (mpc_realref (rop), mpc_realref (rop), GMP_RNDN);
 214               inex_re = 1;
 215             }
 216           else
 217             inex_re = -1; /* +0 is rounded down */
 218           if (mpfr_sgn (mpc_imagref (op)) > 0)
 219             {
 220               mpfr_set_ui (mpc_imagref (rop), 1, GMP_RNDN);
 221               inex_im = 1;
 222             }
 223           else
 224             {
 225               mpfr_set_si (mpc_imagref (rop), -1, GMP_RNDN);
 226               inex_im = -1;
 227             }
 228           inex = MPC_INEX(inex_re, inex_im);
 229           goto end;
 230         }
 231
 232       exr = mpfr_get_exp (mpc_realref (x));
 233       eyr = mpfr_get_exp (mpc_realref (y));
 234       eyi = mpfr_get_exp (mpc_imagref (y));
 235
 236       /* some parts of the quotient may be exact */
 237       inex = mpc_div (x, x, y, MPC_RNDZZ);
 238       /* OP is no pure real nor pure imaginary, so in theory the real and
 239          imaginary parts of its tangent cannot be null. However due to
 240          rouding errors this might happen. Consider for example
 241          tan(1+14*I) = 1.26e-10 + 1.00*I. For small precision sin(op) and
 242          cos(op) differ only by a factor I, thus after mpc_div x = I and
 243          its real part is zero. */
 244       if (mpfr_zero_p (mpc_realref (x)) || mpfr_zero_p (mpc_imagref (x)))
 245         {
 246           err = prec; /* double precision */
 247           continue;
 248         }
 249       if (MPC_INEX_RE (inex))
 250          MPFR_ADD_ONE_ULP (mpc_realref (x));
 251       if (MPC_INEX_IM (inex))
 252          MPFR_ADD_ONE_ULP (mpc_imagref (x));
 253       MPC_ASSERT (mpfr_zero_p (mpc_realref (x)) == 0);
 254       ezr = mpfr_get_exp (mpc_realref (x));
 255
 256       /* FIXME: compute
 257          k = Exp(Re(x))+Exp(Re(y))-2min{Exp(Re(y)), Exp(Im(y))}-Exp(Re(x/y))
 258          avoiding overflow */
 259       k = exr - ezr + MPC_MAX(-eyr, eyr - 2 * eyi);
 260       err = k < 2 ? 7 : (k == 2 ? 8 : (5 + k));
 261
 262       /* Can the real part be rounded? */
 263       ok = (!mpfr_number_p (mpc_realref (x)))
 264            || mpfr_can_round (mpc_realref(x), prec - err, GMP_RNDN, GMP_RNDZ,
 265                       MPC_PREC_RE(rop) + (MPC_RND_RE(rnd) == GMP_RNDN));
 266
 267       if (ok)
 268         {
 269           /* Can the imaginary part be rounded? */
 270           ok = (!mpfr_number_p (mpc_imagref (x)))
 271                || mpfr_can_round (mpc_imagref(x), prec - 6, GMP_RNDN, GMP_RNDZ,
 272                       MPC_PREC_IM(rop) + (MPC_RND_IM(rnd) == GMP_RNDN));
 273         }
 274     }
 275   while (ok == 0);
 276
 277   inex = mpc_set (rop, x, rnd);
 278
 279  end:
 280   mpc_clear (x);
 281   mpc_clear (y);
 282
 283   return inex;
 284 }