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
mpc_tan(mpc_ptr rop,mpc_srcptr op,mpc_rnd_t rnd)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 }
285