1/* mpc_tan -- tangent of a complex number.
2
3Copyright (C) 2008, 2009, 2010, 2011 INRIA
4
5This file is part of GNU MPC.
6
7GNU MPC is free software; you can redistribute it and/or modify it under
8the terms of the GNU Lesser General Public License as published by the
9Free Software Foundation; either version 3 of the License, or (at your
10option) any later version.
11
12GNU MPC is distributed in the hope that it will be useful, but WITHOUT ANY
13WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
14FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
15more details.
16
17You should have received a copy of the GNU Lesser General Public License
18along 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
25int
26mpc_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}
285