1 /* mpc_tan -- tangent of a complex number.
2
3 Copyright (C) 2008, 2009, 2010, 2011, 2012, 2013, 2015, 2020, 2022 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 /* special case where the imaginary part of tan(op) rounds to -1 or 1:
26 return 1 if |Im(tan(op))| > 1, and -1 if |Im(tan(op))| < 1, return 0
27 if we can't decide.
28 The imaginary part is sinh(2*y)/(cos(2*x) + cosh(2*y)) where op = (x,y).
29 */
30 static int
tan_im_cmp_one(mpc_srcptr op)31 tan_im_cmp_one (mpc_srcptr op)
32 {
33 mpfr_t x, c;
34 int ret = 0;
35 mpfr_exp_t expc;
36
37 mpfr_init2 (x, mpfr_get_prec (mpc_realref (op)));
38 mpfr_mul_2exp (x, mpc_realref (op), 1, MPFR_RNDN);
39 mpfr_init2 (c, 32);
40 mpfr_cos (c, x, MPFR_RNDN);
41 /* if cos(2x) >= 0, then |sinh(2y)/(cos(2x)+cosh(2y))| < 1 */
42 if (mpfr_sgn (c) >= 0)
43 ret = -1; /* |Im(tan(op))| < 1 */
44 else
45 {
46 /* now cos(2x) < 0: |cosh(2y) - sinh(2y)| = exp(-2|y|) */
47 expc = mpfr_get_exp (c);
48 mpfr_abs (c, mpc_imagref (op), MPFR_RNDN);
49 mpfr_mul_si (c, c, -2, MPFR_RNDN);
50 mpfr_exp (c, c, MPFR_RNDN);
51 if (mpfr_zero_p (c) || mpfr_get_exp (c) < expc)
52 ret = 1; /* |Im(tan(op))| > 1 */
53 }
54 mpfr_clear (c);
55 mpfr_clear (x);
56 return ret;
57 }
58
59 /* special case where the real part of tan(op) underflows to 0:
60 return 1 if 0 < Re(tan(op)) < 2^(emin-2),
61 -1 if -2^(emin-2) < Re(tan(op))| < 0, and 0 if we can't decide.
62 The real part is sin(2*x)/(cos(2*x) + cosh(2*y)) where op = (x,y),
63 thus has the sign of sin(2*x).
64 */
65 static int
tan_re_cmp_zero(mpc_srcptr op,mpfr_exp_t emin)66 tan_re_cmp_zero (mpc_srcptr op, mpfr_exp_t emin)
67 {
68 mpfr_t x, s, c;
69 int ret = 0;
70
71 mpfr_init2 (x, mpfr_get_prec (mpc_realref (op)));
72 mpfr_mul_2exp (x, mpc_realref (op), 1, MPFR_RNDN);
73 mpfr_init2 (s, 32);
74 mpfr_init2 (c, 32);
75 mpfr_sin (s, x, MPFR_RNDA);
76 mpfr_mul_2exp (x, mpc_imagref (op), 1, MPFR_RNDN);
77 mpfr_cosh (c, x, MPFR_RNDZ);
78 mpfr_sub_ui (c, c, 1, MPFR_RNDZ);
79 mpfr_div (s, s, c, MPFR_RNDA);
80 if (mpfr_zero_p (s) || mpfr_get_exp (s) <= emin - 2)
81 ret = mpfr_sgn (s);
82 mpfr_clear (s);
83 mpfr_clear (c);
84 mpfr_clear (x);
85 return ret;
86 }
87
88 int
mpc_tan(mpc_ptr rop,mpc_srcptr op,mpc_rnd_t rnd)89 mpc_tan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd)
90 {
91 mpc_t x, y;
92 mpfr_prec_t prec;
93 mpfr_exp_t err;
94 int ok;
95 int inex, inex_re, inex_im;
96 mpfr_exp_t saved_emin, saved_emax;
97
98 /* special values */
99 if (!mpc_fin_p (op))
100 {
101 if (mpfr_nan_p (mpc_realref (op)))
102 {
103 if (mpfr_inf_p (mpc_imagref (op)))
104 /* tan(NaN -i*Inf) = +/-0 -i */
105 /* tan(NaN +i*Inf) = +/-0 +i */
106 {
107 /* exact unless 1 is not in exponent range */
108 inex = mpc_set_si_si (rop, 0,
109 (MPFR_SIGN (mpc_imagref (op)) < 0) ? -1 : +1,
110 rnd);
111 }
112 else
113 /* tan(NaN +i*y) = NaN +i*NaN, when y is finite */
114 /* tan(NaN +i*NaN) = NaN +i*NaN */
115 {
116 mpfr_set_nan (mpc_realref (rop));
117 mpfr_set_nan (mpc_imagref (rop));
118 inex = MPC_INEX (0, 0); /* always exact */
119 }
120 }
121 else if (mpfr_nan_p (mpc_imagref (op)))
122 {
123 if (mpfr_cmp_ui (mpc_realref (op), 0) == 0)
124 /* tan(-0 +i*NaN) = -0 +i*NaN */
125 /* tan(+0 +i*NaN) = +0 +i*NaN */
126 {
127 mpc_set (rop, op, rnd);
128 inex = MPC_INEX (0, 0); /* always exact */
129 }
130 else
131 /* tan(x +i*NaN) = NaN +i*NaN, when x != 0 */
132 {
133 mpfr_set_nan (mpc_realref (rop));
134 mpfr_set_nan (mpc_imagref (rop));
135 inex = MPC_INEX (0, 0); /* always exact */
136 }
137 }
138 else if (mpfr_inf_p (mpc_realref (op)))
139 {
140 if (mpfr_inf_p (mpc_imagref (op)))
141 /* tan(-Inf -i*Inf) = -/+0 -i */
142 /* tan(-Inf +i*Inf) = -/+0 +i */
143 /* tan(+Inf -i*Inf) = +/-0 -i */
144 /* tan(+Inf +i*Inf) = +/-0 +i */
145 {
146 const int sign_re = mpfr_signbit (mpc_realref (op));
147
148 mpfr_set_ui (mpc_realref (rop), 0, MPC_RND_RE (rnd));
149 mpfr_setsign (mpc_realref (rop), mpc_realref (rop), sign_re, MPFR_RNDN);
150
151 /* exact, unless 1 is not in exponent range */
152 inex_im = mpfr_set_si (mpc_imagref (rop),
153 mpfr_signbit (mpc_imagref (op)) ? -1 : +1,
154 MPC_RND_IM (rnd));
155 inex = MPC_INEX (0, inex_im);
156 }
157 else
158 /* tan(-Inf +i*y) = tan(+Inf +i*y) = NaN +i*NaN, when y is
159 finite */
160 {
161 mpfr_set_nan (mpc_realref (rop));
162 mpfr_set_nan (mpc_imagref (rop));
163 inex = MPC_INEX (0, 0); /* always exact */
164 }
165 }
166 else
167 /* tan(x -i*Inf) = +0*sin(x)*cos(x) -i, when x is finite */
168 /* tan(x +i*Inf) = +0*sin(x)*cos(x) +i, when x is finite */
169 {
170 mpfr_t c;
171 mpfr_t s;
172
173 mpfr_init (c);
174 mpfr_init (s);
175
176 mpfr_sin_cos (s, c, mpc_realref (op), MPFR_RNDN);
177 mpfr_set_ui (mpc_realref (rop), 0, MPC_RND_RE (rnd));
178 mpfr_setsign (mpc_realref (rop), mpc_realref (rop),
179 mpfr_signbit (c) != mpfr_signbit (s), MPFR_RNDN);
180 /* exact, unless 1 is not in exponent range */
181 inex_im = mpfr_set_si (mpc_imagref (rop),
182 (mpfr_signbit (mpc_imagref (op)) ? -1 : +1),
183 MPC_RND_IM (rnd));
184 inex = MPC_INEX (0, inex_im);
185
186 mpfr_clear (s);
187 mpfr_clear (c);
188 }
189
190 return inex;
191 }
192
193 if (mpfr_zero_p (mpc_realref (op)))
194 /* tan(-0 -i*y) = -0 +i*tanh(y), when y is finite. */
195 /* tan(+0 +i*y) = +0 +i*tanh(y), when y is finite. */
196 {
197 mpfr_set (mpc_realref (rop), mpc_realref (op), MPC_RND_RE (rnd));
198 inex_im = mpfr_tanh (mpc_imagref (rop), mpc_imagref (op), MPC_RND_IM (rnd));
199
200 return MPC_INEX (0, inex_im);
201 }
202
203 if (mpfr_zero_p (mpc_imagref (op)))
204 /* tan(x -i*0) = tan(x) -i*0, when x is finite. */
205 /* tan(x +i*0) = tan(x) +i*0, when x is finite. */
206 {
207 inex_re = mpfr_tan (mpc_realref (rop), mpc_realref (op), MPC_RND_RE (rnd));
208 mpfr_set (mpc_imagref (rop), mpc_imagref (op), MPC_RND_IM (rnd));
209
210 return MPC_INEX (inex_re, 0);
211 }
212
213 saved_emin = mpfr_get_emin ();
214 saved_emax = mpfr_get_emax ();
215 mpfr_set_emin (mpfr_get_emin_min ());
216 mpfr_set_emax (mpfr_get_emax_max ());
217
218 /* ordinary (non-zero) numbers */
219
220 /* tan(op) = sin(op) / cos(op).
221
222 We use the following algorithm with rounding away from 0 for all
223 operations, and working precision w:
224
225 (1) x = A(sin(op))
226 (2) y = A(cos(op))
227 (3) z = A(x/y)
228
229 the error on Im(z) is at most 81 ulp,
230 the error on Re(z) is at most
231 7 ulp if k < 2,
232 8 ulp if k = 2,
233 else 5+k ulp, where
234 k = Exp(Re(x))+Exp(Re(y))-2min{Exp(Re(y)), Exp(Im(y))}-Exp(Re(x/y))
235 see proof in algorithms.tex.
236 */
237
238 prec = MPC_MAX_PREC(rop);
239
240 mpc_init2 (x, 2);
241 mpc_init2 (y, 2);
242
243 err = 7;
244
245 do
246 {
247 mpfr_exp_t k, exr, eyr, eyi, ezr;
248
249 ok = 0;
250
251 /* FIXME: prevent addition overflow */
252 prec += mpc_ceil_log2 (prec) + err;
253 mpc_set_prec (x, prec);
254 mpc_set_prec (y, prec);
255
256 mpc_sin_cos (x, y, op, MPC_RNDAA, MPC_RNDAA);
257
258 if ( mpfr_inf_p (mpc_realref (x)) || mpfr_inf_p (mpc_imagref (x))
259 || mpfr_inf_p (mpc_realref (y)) || mpfr_inf_p (mpc_imagref (y))) {
260 /* If the real or imaginary part of x is infinite, it means that
261 Im(op) was large, in which case the result is
262 sign(tan(Re(op)))*0 + sign(Im(op))*I,
263 where sign(tan(Re(op))) = sign(Re(x))*sign(Re(y)). */
264 mpfr_set_ui (mpc_realref (rop), 0, MPFR_RNDN);
265 if (mpfr_sgn (mpc_realref (x)) * mpfr_sgn (mpc_realref (y)) < 0)
266 {
267 mpfr_neg (mpc_realref (rop), mpc_realref (rop), MPFR_RNDN);
268 inex_re = 1;
269 }
270 else
271 inex_re = -1; /* +0 is rounded down */
272 if (mpfr_sgn (mpc_imagref (op)) > 0)
273 {
274 mpfr_set_ui (mpc_imagref (rop), 1, MPFR_RNDN);
275 inex_im = 1;
276 }
277 else
278 {
279 mpfr_set_si (mpc_imagref (rop), -1, MPFR_RNDN);
280 inex_im = -1;
281 }
282 /* if rounding is toward zero, fix the imaginary part */
283 if (MPC_IS_LIKE_RNDZ(MPC_RND_IM(rnd), MPFR_SIGNBIT(mpc_imagref (rop))))
284 {
285 mpfr_nexttoward (mpc_imagref (rop), mpc_realref (rop));
286 inex_im = -inex_im;
287 }
288 if (mpfr_zero_p (mpc_realref (rop)))
289 inex_re = mpc_fix_zero (mpc_realref (rop), MPC_RND_RE(rnd));
290 if (mpfr_zero_p (mpc_imagref (rop)))
291 inex_im = mpc_fix_zero (mpc_imagref (rop), MPC_RND_IM(rnd));
292 inex = MPC_INEX(inex_re, inex_im);
293 goto end;
294 }
295
296 exr = mpfr_get_exp (mpc_realref (x));
297 eyr = mpfr_get_exp (mpc_realref (y));
298 eyi = mpfr_get_exp (mpc_imagref (y));
299
300 /* some parts of the quotient may be exact */
301 inex = mpc_div (x, x, y, MPC_RNDZZ);
302 /* OP is no pure real nor pure imaginary, so in theory the real and
303 imaginary parts of its tangent cannot be null. However due to
304 rounding errors this might happen. Consider for example
305 tan(1+14*I) = 1.26e-10 + 1.00*I. For small precision sin(op) and
306 cos(op) differ only by a factor I, thus after mpc_div x = I and
307 its real part is zero. */
308 if (mpfr_zero_p (mpc_realref (x)))
309 {
310 /* since we use an extended exponent range, if real(x) is zero,
311 this means the real part underflows, and we assume we can round */
312 ok = tan_re_cmp_zero (op, saved_emin);
313 if (ok > 0)
314 MPFR_ADD_ONE_ULP (mpc_realref (x));
315 else
316 MPFR_SUB_ONE_ULP (mpc_realref (x));
317 }
318 else
319 {
320 if (MPC_INEX_RE (inex))
321 MPFR_ADD_ONE_ULP (mpc_realref (x));
322 MPC_ASSERT (mpfr_zero_p (mpc_realref (x)) == 0);
323 ezr = mpfr_get_exp (mpc_realref (x));
324
325 /* FIXME: compute
326 k = Exp(Re(x))+Exp(Re(y))-2min{Exp(Re(y)), Exp(Im(y))}-Exp(Re(x/y))
327 avoiding overflow */
328 k = exr - ezr + MPC_MAX(-eyr, eyr - 2 * eyi);
329 err = k < 2 ? 7 : (k == 2 ? 8 : (5 + k));
330
331 /* Can the real part be rounded? */
332 ok = (!mpfr_number_p (mpc_realref (x)))
333 || mpfr_can_round (mpc_realref(x), prec - err, MPFR_RNDN, MPFR_RNDZ,
334 MPC_PREC_RE(rop) + (MPC_RND_RE(rnd) == MPFR_RNDN));
335 }
336
337 if (ok)
338 {
339 if (MPC_INEX_IM (inex))
340 MPFR_ADD_ONE_ULP (mpc_imagref (x));
341
342 /* Can the imaginary part be rounded? */
343 ok = (!mpfr_number_p (mpc_imagref (x)))
344 || mpfr_can_round (mpc_imagref(x), prec - 6, MPFR_RNDN, MPFR_RNDZ,
345 MPC_PREC_IM(rop) + (MPC_RND_IM(rnd) == MPFR_RNDN));
346
347 /* Special case when Im(x) = +/- 1:
348 tan z = [sin(2x)+i*sinh(2y)] / [cos(2x) + cosh(2y)]
349 (formula 4.3.57 of Abramowitz and Stegun) thus for y large
350 in absolute value the imaginary part is near -1 or +1.
351 More precisely cos(2x) + cosh(2y) = cosh(2y) + t with |t| <= 1,
352 thus since cosh(2y) >= exp|2y|/2, then the imaginary part is:
353 tanh(2y) * 1/(1+u) where u = |cos(2x)/cosh(2y)| <= 2/exp|2y|
354 thus |im(z) - tanh(2y)| <= 2/exp|2y| * tanh(2y).
355 Since |tanh(2y)| = (1-exp(-4|y|))/(1+exp(-4|y|)),
356 we have 1-|tanh(2y)| < 2*exp(-4|y|).
357 Thus |im(z)-1| < 2/exp|2y| + 2/exp|4y| < 4/exp|2y| < 4/2^|2y|.
358 If 2^EXP(y) >= p+2, then im(z) rounds to -1 or 1. */
359 if (ok == 0 && (mpfr_cmp_ui (mpc_imagref(x), 1) == 0 ||
360 mpfr_cmp_si (mpc_imagref(x), -1) == 0) &&
361 mpfr_get_exp (mpc_imagref(op)) >= 0 &&
362 ((size_t) mpfr_get_exp (mpc_imagref(op)) >= 8 * sizeof (mpfr_prec_t) ||
363 ((mpfr_prec_t) 1) << mpfr_get_exp (mpc_imagref(op)) >= mpfr_get_prec (mpc_imagref (rop)) + 2))
364 {
365 /* subtract one ulp, so that we get the correct inexact flag */
366 ok = tan_im_cmp_one (op);
367 if (ok < 0)
368 MPFR_SUB_ONE_ULP (mpc_imagref(x));
369 else if (ok > 0)
370 MPFR_ADD_ONE_ULP (mpc_imagref(x));
371 }
372 }
373
374 if (ok == 0)
375 prec += prec / 2;
376 }
377 while (ok == 0);
378
379 inex = mpc_set (rop, x, rnd);
380
381 end:
382 mpc_clear (x);
383 mpc_clear (y);
384
385 /* restore the exponent range, and check the range of results */
386 mpfr_set_emin (saved_emin);
387 mpfr_set_emax (saved_emax);
388 inex_re = mpfr_check_range (mpc_realref (rop), MPC_INEX_RE(inex),
389 MPC_RND_RE (rnd));
390 inex_im = mpfr_check_range (mpc_imagref (rop), MPC_INEX_IM(inex),
391 MPC_RND_IM (rnd));
392
393 return MPC_INEX(inex_re, inex_im);
394 }
395