1 /* mpc_sin_cos -- combined sine and cosine of a complex number.
2
3 Copyright (C) 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 "mpc-impl.h"
22
23 static int
mpc_sin_cos_nonfinite(mpc_ptr rop_sin,mpc_ptr rop_cos,mpc_srcptr op,mpc_rnd_t rnd_sin,mpc_rnd_t rnd_cos)24 mpc_sin_cos_nonfinite (mpc_ptr rop_sin, mpc_ptr rop_cos, mpc_srcptr op,
25 mpc_rnd_t rnd_sin, mpc_rnd_t rnd_cos)
26 /* assumes that op (that is, its real or imaginary part) is not finite */
27 {
28 int overlap;
29 mpc_t op_loc;
30
31 overlap = (rop_sin == op || rop_cos == op);
32 if (overlap) {
33 mpc_init3 (op_loc, MPC_PREC_RE (op), MPC_PREC_IM (op));
34 mpc_set (op_loc, op, MPC_RNDNN);
35 }
36 else
37 op_loc [0] = op [0];
38
39 if (rop_sin != NULL) {
40 if (mpfr_nan_p (mpc_realref (op_loc)) || mpfr_nan_p (mpc_imagref (op_loc))) {
41 mpc_set (rop_sin, op_loc, rnd_sin);
42 if (mpfr_nan_p (mpc_imagref (op_loc))) {
43 /* sin(x +i*NaN) = NaN +i*NaN, except for x=0 */
44 /* sin(-0 +i*NaN) = -0 +i*NaN */
45 /* sin(+0 +i*NaN) = +0 +i*NaN */
46 if (!mpfr_zero_p (mpc_realref (op_loc)))
47 mpfr_set_nan (mpc_realref (rop_sin));
48 }
49 else /* op = NaN + i*y */
50 if (!mpfr_inf_p (mpc_imagref (op_loc)) && !mpfr_zero_p (mpc_imagref (op_loc)))
51 /* sin(NaN -i*Inf) = NaN -i*Inf */
52 /* sin(NaN -i*0) = NaN -i*0 */
53 /* sin(NaN +i*0) = NaN +i*0 */
54 /* sin(NaN +i*Inf) = NaN +i*Inf */
55 /* sin(NaN +i*y) = NaN +i*NaN, when 0<|y|<Inf */
56 mpfr_set_nan (mpc_imagref (rop_sin));
57 }
58 else if (mpfr_inf_p (mpc_realref (op_loc))) {
59 mpfr_set_nan (mpc_realref (rop_sin));
60
61 if (!mpfr_inf_p (mpc_imagref (op_loc)) && !mpfr_zero_p (mpc_imagref (op_loc)))
62 /* sin(+/-Inf +i*y) = NaN +i*NaN, when 0<|y|<Inf */
63 mpfr_set_nan (mpc_imagref (rop_sin));
64 else
65 /* sin(+/-Inf -i*Inf) = NaN -i*Inf */
66 /* sin(+/-Inf +i*Inf) = NaN +i*Inf */
67 /* sin(+/-Inf -i*0) = NaN -i*0 */
68 /* sin(+/-Inf +i*0) = NaN +i*0 */
69 mpfr_set (mpc_imagref (rop_sin), mpc_imagref (op_loc), MPC_RND_IM (rnd_sin));
70 }
71 else if (mpfr_zero_p (mpc_realref (op_loc))) {
72 /* sin(-0 -i*Inf) = -0 -i*Inf */
73 /* sin(+0 -i*Inf) = +0 -i*Inf */
74 /* sin(-0 +i*Inf) = -0 +i*Inf */
75 /* sin(+0 +i*Inf) = +0 +i*Inf */
76 mpc_set (rop_sin, op_loc, rnd_sin);
77 }
78 else {
79 /* sin(x -i*Inf) = +Inf*(sin(x) -i*cos(x)) */
80 /* sin(x +i*Inf) = +Inf*(sin(x) +i*cos(x)) */
81 mpfr_t s, c;
82 mpfr_init2 (s, 2);
83 mpfr_init2 (c, 2);
84 mpfr_sin_cos (s, c, mpc_realref (op_loc), GMP_RNDZ);
85 mpfr_set_inf (mpc_realref (rop_sin), MPFR_SIGN (s));
86 mpfr_set_inf (mpc_imagref (rop_sin), MPFR_SIGN (c)*MPFR_SIGN (mpc_imagref (op_loc)));
87 mpfr_clear (s);
88 mpfr_clear (c);
89 }
90 }
91
92 if (rop_cos != NULL) {
93 if (mpfr_nan_p (mpc_realref (op_loc))) {
94 /* cos(NaN + i * NaN) = NaN + i * NaN */
95 /* cos(NaN - i * Inf) = +Inf + i * NaN */
96 /* cos(NaN + i * Inf) = +Inf + i * NaN */
97 /* cos(NaN - i * 0) = NaN - i * 0 */
98 /* cos(NaN + i * 0) = NaN + i * 0 */
99 /* cos(NaN + i * y) = NaN + i * NaN, when y != 0 */
100 if (mpfr_inf_p (mpc_imagref (op_loc)))
101 mpfr_set_inf (mpc_realref (rop_cos), +1);
102 else
103 mpfr_set_nan (mpc_realref (rop_cos));
104
105 if (mpfr_zero_p (mpc_imagref (op_loc)))
106 mpfr_set (mpc_imagref (rop_cos), mpc_imagref (op_loc), MPC_RND_IM (rnd_cos));
107 else
108 mpfr_set_nan (mpc_imagref (rop_cos));
109 }
110 else if (mpfr_nan_p (mpc_imagref (op_loc))) {
111 /* cos(-Inf + i * NaN) = NaN + i * NaN */
112 /* cos(+Inf + i * NaN) = NaN + i * NaN */
113 /* cos(-0 + i * NaN) = NaN - i * 0 */
114 /* cos(+0 + i * NaN) = NaN + i * 0 */
115 /* cos(x + i * NaN) = NaN + i * NaN, when x != 0 */
116 if (mpfr_zero_p (mpc_realref (op_loc)))
117 mpfr_set (mpc_imagref (rop_cos), mpc_realref (op_loc), MPC_RND_IM (rnd_cos));
118 else
119 mpfr_set_nan (mpc_imagref (rop_cos));
120
121 mpfr_set_nan (mpc_realref (rop_cos));
122 }
123 else if (mpfr_inf_p (mpc_realref (op_loc))) {
124 /* cos(-Inf -i*Inf) = cos(+Inf +i*Inf) = -Inf +i*NaN */
125 /* cos(-Inf +i*Inf) = cos(+Inf -i*Inf) = +Inf +i*NaN */
126 /* cos(-Inf -i*0) = cos(+Inf +i*0) = NaN -i*0 */
127 /* cos(-Inf +i*0) = cos(+Inf -i*0) = NaN +i*0 */
128 /* cos(-Inf +i*y) = cos(+Inf +i*y) = NaN +i*NaN, when y != 0 */
129
130 const int same_sign =
131 mpfr_signbit (mpc_realref (op_loc)) == mpfr_signbit (mpc_imagref (op_loc));
132
133 if (mpfr_inf_p (mpc_imagref (op_loc)))
134 mpfr_set_inf (mpc_realref (rop_cos), (same_sign ? -1 : +1));
135 else
136 mpfr_set_nan (mpc_realref (rop_cos));
137
138 if (mpfr_zero_p (mpc_imagref (op_loc)))
139 mpfr_setsign (mpc_imagref (rop_cos), mpc_imagref (op_loc), same_sign,
140 MPC_RND_IM(rnd_cos));
141 else
142 mpfr_set_nan (mpc_imagref (rop_cos));
143 }
144 else if (mpfr_zero_p (mpc_realref (op_loc))) {
145 /* cos(-0 -i*Inf) = cos(+0 +i*Inf) = +Inf -i*0 */
146 /* cos(-0 +i*Inf) = cos(+0 -i*Inf) = +Inf +i*0 */
147 const int same_sign =
148 mpfr_signbit (mpc_realref (op_loc)) == mpfr_signbit (mpc_imagref (op_loc));
149
150 mpfr_setsign (mpc_imagref (rop_cos), mpc_realref (op_loc), same_sign,
151 MPC_RND_IM (rnd_cos));
152 mpfr_set_inf (mpc_realref (rop_cos), +1);
153 }
154 else {
155 /* cos(x -i*Inf) = +Inf*cos(x) +i*Inf*sin(x), when x != 0 */
156 /* cos(x +i*Inf) = +Inf*cos(x) -i*Inf*sin(x), when x != 0 */
157 mpfr_t s, c;
158 mpfr_init2 (c, 2);
159 mpfr_init2 (s, 2);
160 mpfr_sin_cos (s, c, mpc_realref (op_loc), GMP_RNDN);
161 mpfr_set_inf (mpc_realref (rop_cos), mpfr_sgn (c));
162 mpfr_set_inf (mpc_imagref (rop_cos),
163 (mpfr_sgn (mpc_imagref (op_loc)) == mpfr_sgn (s) ? -1 : +1));
164 mpfr_clear (s);
165 mpfr_clear (c);
166 }
167 }
168
169 if (overlap)
170 mpc_clear (op_loc);
171
172 return MPC_INEX12 (MPC_INEX (0,0), MPC_INEX (0,0));
173 /* everything is exact */
174 }
175
176
177 static int
mpc_sin_cos_real(mpc_ptr rop_sin,mpc_ptr rop_cos,mpc_srcptr op,mpc_rnd_t rnd_sin,mpc_rnd_t rnd_cos)178 mpc_sin_cos_real (mpc_ptr rop_sin, mpc_ptr rop_cos, mpc_srcptr op,
179 mpc_rnd_t rnd_sin, mpc_rnd_t rnd_cos)
180 /* assumes that op is real */
181 {
182 int inex_sin_re = 0, inex_cos_re = 0;
183 /* Until further notice, assume computations exact; in particular,
184 by definition, for not computed values. */
185 mpfr_t s, c;
186 int inex_s, inex_c;
187 int sign_im = mpfr_signbit (mpc_imagref (op));
188
189 /* sin(x +-0*i) = sin(x) +-0*i*sign(cos(x)) */
190 /* cos(x +-i*0) = cos(x) -+i*0*sign(sin(x)) */
191 if (rop_sin != 0)
192 mpfr_init2 (s, MPC_PREC_RE (rop_sin));
193 else
194 mpfr_init2 (s, 2); /* We need only the sign. */
195 if (rop_cos != NULL)
196 mpfr_init2 (c, MPC_PREC_RE (rop_cos));
197 else
198 mpfr_init2 (c, 2);
199 inex_s = mpfr_sin (s, mpc_realref (op), MPC_RND_RE (rnd_sin));
200 inex_c = mpfr_cos (c, mpc_realref (op), MPC_RND_RE (rnd_cos));
201 /* We cannot use mpfr_sin_cos since we may need two distinct rounding
202 modes and the exact return values. If we need only the sign, an
203 arbitrary rounding mode will work. */
204
205 if (rop_sin != NULL) {
206 mpfr_set (mpc_realref (rop_sin), s, GMP_RNDN); /* exact */
207 inex_sin_re = inex_s;
208 mpfr_set_zero (mpc_imagref (rop_sin),
209 ( ( sign_im && !mpfr_signbit(c))
210 || (!sign_im && mpfr_signbit(c)) ? -1 : 1));
211 }
212
213 if (rop_cos != NULL) {
214 mpfr_set (mpc_realref (rop_cos), c, GMP_RNDN); /* exact */
215 inex_cos_re = inex_c;
216 mpfr_set_zero (mpc_imagref (rop_cos),
217 ( ( sign_im && mpfr_signbit(s))
218 || (!sign_im && !mpfr_signbit(s)) ? -1 : 1));
219 }
220
221 mpfr_clear (s);
222 mpfr_clear (c);
223
224 return MPC_INEX12 (MPC_INEX (inex_sin_re, 0), MPC_INEX (inex_cos_re, 0));
225 }
226
227
228 static int
mpc_sin_cos_imag(mpc_ptr rop_sin,mpc_ptr rop_cos,mpc_srcptr op,mpc_rnd_t rnd_sin,mpc_rnd_t rnd_cos)229 mpc_sin_cos_imag (mpc_ptr rop_sin, mpc_ptr rop_cos, mpc_srcptr op,
230 mpc_rnd_t rnd_sin, mpc_rnd_t rnd_cos)
231 /* assumes that op is purely imaginary, but not zero */
232 {
233 int inex_sin_im = 0, inex_cos_re = 0;
234 /* assume exact if not computed */
235 int overlap;
236 mpc_t op_loc;
237
238 overlap = (rop_sin == op || rop_cos == op);
239 if (overlap) {
240 mpc_init3 (op_loc, MPC_PREC_RE (op), MPC_PREC_IM (op));
241 mpc_set (op_loc, op, MPC_RNDNN);
242 }
243 else
244 op_loc [0] = op [0];
245
246 if (rop_sin != NULL) {
247 /* sin(+-O +i*y) = +-0 +i*sinh(y) */
248 mpfr_set (mpc_realref(rop_sin), mpc_realref(op_loc), GMP_RNDN);
249 inex_sin_im = mpfr_sinh (mpc_imagref(rop_sin), mpc_imagref(op_loc), MPC_RND_IM(rnd_sin));
250 }
251
252 if (rop_cos != NULL) {
253 /* cos(-0 - i * y) = cos(+0 + i * y) = cosh(y) - i * 0,
254 cos(-0 + i * y) = cos(+0 - i * y) = cosh(y) + i * 0,
255 where y > 0 */
256 inex_cos_re = mpfr_cosh (mpc_realref (rop_cos), mpc_imagref (op_loc), MPC_RND_RE (rnd_cos));
257
258 mpfr_set_ui (mpc_imagref (rop_cos), 0ul, MPC_RND_IM (rnd_cos));
259 if (mpfr_signbit (mpc_realref (op_loc)) == mpfr_signbit (mpc_imagref (op_loc)))
260 MPFR_CHANGE_SIGN (mpc_imagref (rop_cos));
261 }
262
263 if (overlap)
264 mpc_clear (op_loc);
265
266 return MPC_INEX12 (MPC_INEX (0, inex_sin_im), MPC_INEX (inex_cos_re, 0));
267 }
268
269
270 int
mpc_sin_cos(mpc_ptr rop_sin,mpc_ptr rop_cos,mpc_srcptr op,mpc_rnd_t rnd_sin,mpc_rnd_t rnd_cos)271 mpc_sin_cos (mpc_ptr rop_sin, mpc_ptr rop_cos, mpc_srcptr op,
272 mpc_rnd_t rnd_sin, mpc_rnd_t rnd_cos)
273 /* Feature not documented in the texinfo file: One of rop_sin or
274 rop_cos may be NULL, in which case it is not computed, and the
275 corresponding ternary inexact value is set to 0 (exact). */
276 {
277 if (!mpc_fin_p (op))
278 return mpc_sin_cos_nonfinite (rop_sin, rop_cos, op, rnd_sin, rnd_cos);
279 else if (mpfr_zero_p (mpc_imagref (op)))
280 return mpc_sin_cos_real (rop_sin, rop_cos, op, rnd_sin, rnd_cos);
281 else if (mpfr_zero_p (mpc_realref (op)))
282 return mpc_sin_cos_imag (rop_sin, rop_cos, op, rnd_sin, rnd_cos);
283 else {
284 /* let op = a + i*b, then sin(op) = sin(a)*cosh(b) + i*cos(a)*sinh(b)
285 and cos(op) = cos(a)*cosh(b) - i*sin(a)*sinh(b).
286
287 For Re(sin(op)) (and analogously, the other parts), we use the
288 following algorithm, with rounding to nearest for all operations
289 and working precision w:
290
291 (1) x = o(sin(a))
292 (2) y = o(cosh(b))
293 (3) r = o(x*y)
294 then the error on r is at most 4 ulps, since we can write
295 r = sin(a)*cosh(b)*(1+t)^3 with |t| <= 2^(-w),
296 thus for w >= 2, r = sin(a)*cosh(b)*(1+4*t) with |t| <= 2^(-w),
297 thus the relative error is bounded by 4*2^(-w) <= 4*ulp(r).
298 */
299 mpfr_t s, c, sh, ch, sch, csh;
300 mpfr_prec_t prec;
301 int ok;
302 int inex_re, inex_im, inex_sin, inex_cos;
303
304 prec = 2;
305 if (rop_sin != NULL)
306 prec = MPC_MAX (prec, MPC_MAX_PREC (rop_sin));
307 if (rop_cos != NULL)
308 prec = MPC_MAX (prec, MPC_MAX_PREC (rop_cos));
309
310 mpfr_init2 (s, 2);
311 mpfr_init2 (c, 2);
312 mpfr_init2 (sh, 2);
313 mpfr_init2 (ch, 2);
314 mpfr_init2 (sch, 2);
315 mpfr_init2 (csh, 2);
316
317 do {
318 ok = 1;
319 prec += mpc_ceil_log2 (prec) + 5;
320
321 mpfr_set_prec (s, prec);
322 mpfr_set_prec (c, prec);
323 mpfr_set_prec (sh, prec);
324 mpfr_set_prec (ch, prec);
325 mpfr_set_prec (sch, prec);
326 mpfr_set_prec (csh, prec);
327
328 mpfr_sin_cos (s, c, mpc_realref(op), GMP_RNDN);
329 mpfr_sinh_cosh (sh, ch, mpc_imagref(op), GMP_RNDN);
330
331 if (rop_sin != NULL) {
332 /* real part of sine */
333 mpfr_mul (sch, s, ch, GMP_RNDN);
334 ok = (!mpfr_number_p (sch))
335 || mpfr_can_round (sch, prec - 2, GMP_RNDN, GMP_RNDZ,
336 MPC_PREC_RE (rop_sin)
337 + (MPC_RND_RE (rnd_sin) == GMP_RNDN));
338
339 if (ok) {
340 /* imaginary part of sine */
341 mpfr_mul (csh, c, sh, GMP_RNDN);
342 ok = (!mpfr_number_p (csh))
343 || mpfr_can_round (csh, prec - 2, GMP_RNDN, GMP_RNDZ,
344 MPC_PREC_IM (rop_sin)
345 + (MPC_RND_IM (rnd_sin) == GMP_RNDN));
346 }
347 }
348
349 if (rop_cos != NULL && ok) {
350 /* real part of cosine */
351 mpfr_mul (c, c, ch, GMP_RNDN);
352 ok = (!mpfr_number_p (c))
353 || mpfr_can_round (c, prec - 2, GMP_RNDN, GMP_RNDZ,
354 MPC_PREC_RE (rop_cos)
355 + (MPC_RND_RE (rnd_cos) == GMP_RNDN));
356
357 if (ok) {
358 /* imaginary part of cosine */
359 mpfr_mul (s, s, sh, GMP_RNDN);
360 mpfr_neg (s, s, GMP_RNDN);
361 ok = (!mpfr_number_p (s))
362 || mpfr_can_round (s, prec - 2, GMP_RNDN, GMP_RNDZ,
363 MPC_PREC_IM (rop_cos)
364 + (MPC_RND_IM (rnd_cos) == GMP_RNDN));
365 }
366 }
367 } while (ok == 0);
368
369 if (rop_sin != NULL) {
370 inex_re = mpfr_set (mpc_realref (rop_sin), sch, MPC_RND_RE (rnd_sin));
371 if (mpfr_inf_p (sch))
372 inex_re = mpfr_sgn (sch);
373 inex_im = mpfr_set (mpc_imagref (rop_sin), csh, MPC_RND_IM (rnd_sin));
374 if (mpfr_inf_p (csh))
375 inex_im = mpfr_sgn (csh);
376 inex_sin = MPC_INEX (inex_re, inex_im);
377 }
378 else
379 inex_sin = MPC_INEX (0,0); /* return exact if not computed */
380
381 if (rop_cos != NULL) {
382 inex_re = mpfr_set (mpc_realref (rop_cos), c, MPC_RND_RE (rnd_cos));
383 if (mpfr_inf_p (c))
384 inex_re = mpfr_sgn (c);
385 inex_im = mpfr_set (mpc_imagref (rop_cos), s, MPC_RND_IM (rnd_cos));
386 if (mpfr_inf_p (s))
387 inex_im = mpfr_sgn (s);
388 inex_cos = MPC_INEX (inex_re, inex_im);
389 }
390 else
391 inex_cos = MPC_INEX (0,0); /* return exact if not computed */
392
393 mpfr_clear (s);
394 mpfr_clear (c);
395 mpfr_clear (sh);
396 mpfr_clear (ch);
397 mpfr_clear (sch);
398 mpfr_clear (csh);
399
400 return (MPC_INEX12 (inex_sin, inex_cos));
401 }
402 }
403