1 /* mpc_sin_cos -- combined sine and cosine of a complex number.
2
3 Copyright (C) 2010, 2011, 2012, 2020 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>
22 #include "mpc-impl.h"
23
24 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)25 mpc_sin_cos_nonfinite (mpc_ptr rop_sin, mpc_ptr rop_cos, mpc_srcptr op,
26 mpc_rnd_t rnd_sin, mpc_rnd_t rnd_cos)
27 /* assumes that op (that is, its real or imaginary part) is not finite */
28 {
29 int overlap;
30 mpc_t op_loc;
31
32 overlap = (rop_sin == op || rop_cos == op);
33 if (overlap) {
34 mpc_init3 (op_loc, MPC_PREC_RE (op), MPC_PREC_IM (op));
35 mpc_set (op_loc, op, MPC_RNDNN);
36 }
37 else
38 op_loc [0] = op [0];
39
40 if (rop_sin != NULL) {
41 if (mpfr_nan_p (mpc_realref (op_loc)) || mpfr_nan_p (mpc_imagref (op_loc))) {
42 mpc_set (rop_sin, op_loc, rnd_sin);
43 if (mpfr_nan_p (mpc_imagref (op_loc))) {
44 /* sin(x +i*NaN) = NaN +i*NaN, except for x=0 */
45 /* sin(-0 +i*NaN) = -0 +i*NaN */
46 /* sin(+0 +i*NaN) = +0 +i*NaN */
47 if (!mpfr_zero_p (mpc_realref (op_loc)))
48 mpfr_set_nan (mpc_realref (rop_sin));
49 }
50 else /* op = NaN + i*y */
51 if (!mpfr_inf_p (mpc_imagref (op_loc)) && !mpfr_zero_p (mpc_imagref (op_loc)))
52 /* sin(NaN -i*Inf) = NaN -i*Inf */
53 /* sin(NaN -i*0) = NaN -i*0 */
54 /* sin(NaN +i*0) = NaN +i*0 */
55 /* sin(NaN +i*Inf) = NaN +i*Inf */
56 /* sin(NaN +i*y) = NaN +i*NaN, when 0<|y|<Inf */
57 mpfr_set_nan (mpc_imagref (rop_sin));
58 }
59 else if (mpfr_inf_p (mpc_realref (op_loc))) {
60 mpfr_set_nan (mpc_realref (rop_sin));
61
62 if (!mpfr_inf_p (mpc_imagref (op_loc)) && !mpfr_zero_p (mpc_imagref (op_loc)))
63 /* sin(+/-Inf +i*y) = NaN +i*NaN, when 0<|y|<Inf */
64 mpfr_set_nan (mpc_imagref (rop_sin));
65 else
66 /* sin(+/-Inf -i*Inf) = NaN -i*Inf */
67 /* sin(+/-Inf +i*Inf) = NaN +i*Inf */
68 /* sin(+/-Inf -i*0) = NaN -i*0 */
69 /* sin(+/-Inf +i*0) = NaN +i*0 */
70 mpfr_set (mpc_imagref (rop_sin), mpc_imagref (op_loc), MPC_RND_IM (rnd_sin));
71 }
72 else if (mpfr_zero_p (mpc_realref (op_loc))) {
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 /* sin(+0 +i*Inf) = +0 +i*Inf */
77 mpc_set (rop_sin, op_loc, rnd_sin);
78 }
79 else {
80 /* sin(x -i*Inf) = +Inf*(sin(x) -i*cos(x)) */
81 /* sin(x +i*Inf) = +Inf*(sin(x) +i*cos(x)) */
82 mpfr_t s, c;
83 mpfr_init2 (s, 2);
84 mpfr_init2 (c, 2);
85 mpfr_sin_cos (s, c, mpc_realref (op_loc), MPFR_RNDZ);
86 mpfr_set_inf (mpc_realref (rop_sin), MPFR_SIGN (s));
87 mpfr_set_inf (mpc_imagref (rop_sin), MPFR_SIGN (c)*MPFR_SIGN (mpc_imagref (op_loc)));
88 mpfr_clear (s);
89 mpfr_clear (c);
90 }
91 }
92
93 if (rop_cos != NULL) {
94 if (mpfr_nan_p (mpc_realref (op_loc))) {
95 /* cos(NaN + i * NaN) = NaN + i * NaN */
96 /* cos(NaN - i * Inf) = +Inf + i * NaN */
97 /* cos(NaN + i * Inf) = +Inf + i * NaN */
98 /* cos(NaN - i * 0) = NaN - i * 0 */
99 /* cos(NaN + i * 0) = NaN + i * 0 */
100 /* cos(NaN + i * y) = NaN + i * NaN, when y != 0 */
101 if (mpfr_inf_p (mpc_imagref (op_loc)))
102 mpfr_set_inf (mpc_realref (rop_cos), +1);
103 else
104 mpfr_set_nan (mpc_realref (rop_cos));
105
106 if (mpfr_zero_p (mpc_imagref (op_loc)))
107 mpfr_set (mpc_imagref (rop_cos), mpc_imagref (op_loc), MPC_RND_IM (rnd_cos));
108 else
109 mpfr_set_nan (mpc_imagref (rop_cos));
110 }
111 else if (mpfr_nan_p (mpc_imagref (op_loc))) {
112 /* cos(-Inf + i * NaN) = NaN + i * NaN */
113 /* cos(+Inf + i * NaN) = NaN + i * NaN */
114 /* cos(-0 + i * NaN) = NaN - i * 0 */
115 /* cos(+0 + i * NaN) = NaN + i * 0 */
116 /* cos(x + i * NaN) = NaN + i * NaN, when x != 0 */
117 if (mpfr_zero_p (mpc_realref (op_loc)))
118 mpfr_set (mpc_imagref (rop_cos), mpc_realref (op_loc), MPC_RND_IM (rnd_cos));
119 else
120 mpfr_set_nan (mpc_imagref (rop_cos));
121
122 mpfr_set_nan (mpc_realref (rop_cos));
123 }
124 else if (mpfr_inf_p (mpc_realref (op_loc))) {
125 /* cos(-Inf -i*Inf) = cos(+Inf +i*Inf) = -Inf +i*NaN */
126 /* cos(-Inf +i*Inf) = cos(+Inf -i*Inf) = +Inf +i*NaN */
127 /* cos(-Inf -i*0) = cos(+Inf +i*0) = NaN -i*0 */
128 /* cos(-Inf +i*0) = cos(+Inf -i*0) = NaN +i*0 */
129 /* cos(-Inf +i*y) = cos(+Inf +i*y) = NaN +i*NaN, when y != 0 */
130
131 const int same_sign =
132 mpfr_signbit (mpc_realref (op_loc)) == mpfr_signbit (mpc_imagref (op_loc));
133
134 if (mpfr_inf_p (mpc_imagref (op_loc)))
135 mpfr_set_inf (mpc_realref (rop_cos), (same_sign ? -1 : +1));
136 else
137 mpfr_set_nan (mpc_realref (rop_cos));
138
139 if (mpfr_zero_p (mpc_imagref (op_loc)))
140 mpfr_setsign (mpc_imagref (rop_cos), mpc_imagref (op_loc), same_sign,
141 MPC_RND_IM(rnd_cos));
142 else
143 mpfr_set_nan (mpc_imagref (rop_cos));
144 }
145 else if (mpfr_zero_p (mpc_realref (op_loc))) {
146 /* cos(-0 -i*Inf) = cos(+0 +i*Inf) = +Inf -i*0 */
147 /* cos(-0 +i*Inf) = cos(+0 -i*Inf) = +Inf +i*0 */
148 const int same_sign =
149 mpfr_signbit (mpc_realref (op_loc)) == mpfr_signbit (mpc_imagref (op_loc));
150
151 mpfr_setsign (mpc_imagref (rop_cos), mpc_realref (op_loc), same_sign,
152 MPC_RND_IM (rnd_cos));
153 mpfr_set_inf (mpc_realref (rop_cos), +1);
154 }
155 else {
156 /* cos(x -i*Inf) = +Inf*cos(x) +i*Inf*sin(x), when x != 0 */
157 /* cos(x +i*Inf) = +Inf*cos(x) -i*Inf*sin(x), when x != 0 */
158 mpfr_t s, c;
159 mpfr_init2 (c, 2);
160 mpfr_init2 (s, 2);
161 mpfr_sin_cos (s, c, mpc_realref (op_loc), MPFR_RNDN);
162 mpfr_set_inf (mpc_realref (rop_cos), mpfr_sgn (c));
163 mpfr_set_inf (mpc_imagref (rop_cos),
164 (mpfr_sgn (mpc_imagref (op_loc)) == mpfr_sgn (s) ? -1 : +1));
165 mpfr_clear (s);
166 mpfr_clear (c);
167 }
168 }
169
170 if (overlap)
171 mpc_clear (op_loc);
172
173 return MPC_INEX12 (MPC_INEX (0,0), MPC_INEX (0,0));
174 /* everything is exact */
175 }
176
177
178 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)179 mpc_sin_cos_real (mpc_ptr rop_sin, mpc_ptr rop_cos, mpc_srcptr op,
180 mpc_rnd_t rnd_sin, mpc_rnd_t rnd_cos)
181 /* assumes that op is real */
182 {
183 int inex_sin_re = 0, inex_cos_re = 0;
184 /* Until further notice, assume computations exact; in particular,
185 by definition, for not computed values. */
186 mpfr_t s, c;
187 int inex_s, inex_c;
188 int sign_im = mpfr_signbit (mpc_imagref (op));
189
190 /* sin(x +-0*i) = sin(x) +-0*i*sign(cos(x)) */
191 /* cos(x +-i*0) = cos(x) -+i*0*sign(sin(x)) */
192 if (rop_sin != 0)
193 mpfr_init2 (s, MPC_PREC_RE (rop_sin));
194 else
195 mpfr_init2 (s, 2); /* We need only the sign. */
196 if (rop_cos != NULL)
197 mpfr_init2 (c, MPC_PREC_RE (rop_cos));
198 else
199 mpfr_init2 (c, 2);
200 inex_s = mpfr_sin (s, mpc_realref (op), MPC_RND_RE (rnd_sin));
201 inex_c = mpfr_cos (c, mpc_realref (op), MPC_RND_RE (rnd_cos));
202 /* We cannot use mpfr_sin_cos since we may need two distinct rounding
203 modes and the exact return values. If we need only the sign, an
204 arbitrary rounding mode will work. */
205
206 if (rop_sin != NULL) {
207 mpfr_set (mpc_realref (rop_sin), s, MPFR_RNDN); /* exact */
208 inex_sin_re = inex_s;
209 mpfr_set_zero (mpc_imagref (rop_sin),
210 ( ( sign_im && !mpfr_signbit(c))
211 || (!sign_im && mpfr_signbit(c)) ? -1 : 1));
212 }
213
214 if (rop_cos != NULL) {
215 mpfr_set (mpc_realref (rop_cos), c, MPFR_RNDN); /* exact */
216 inex_cos_re = inex_c;
217 mpfr_set_zero (mpc_imagref (rop_cos),
218 ( ( sign_im && mpfr_signbit(s))
219 || (!sign_im && !mpfr_signbit(s)) ? -1 : 1));
220 }
221
222 mpfr_clear (s);
223 mpfr_clear (c);
224
225 return MPC_INEX12 (MPC_INEX (inex_sin_re, 0), MPC_INEX (inex_cos_re, 0));
226 }
227
228
229 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)230 mpc_sin_cos_imag (mpc_ptr rop_sin, mpc_ptr rop_cos, mpc_srcptr op,
231 mpc_rnd_t rnd_sin, mpc_rnd_t rnd_cos)
232 /* assumes that op is purely imaginary, but not zero */
233 {
234 int inex_sin_im = 0, inex_cos_re = 0;
235 /* assume exact if not computed */
236 int overlap;
237 mpc_t op_loc;
238
239 overlap = (rop_sin == op || rop_cos == op);
240 if (overlap) {
241 mpc_init3 (op_loc, MPC_PREC_RE (op), MPC_PREC_IM (op));
242 mpc_set (op_loc, op, MPC_RNDNN);
243 }
244 else
245 op_loc [0] = op [0];
246
247 if (rop_sin != NULL) {
248 /* sin(+-O +i*y) = +-0 +i*sinh(y) */
249 mpfr_set (mpc_realref(rop_sin), mpc_realref(op_loc), MPFR_RNDN);
250 inex_sin_im = mpfr_sinh (mpc_imagref(rop_sin), mpc_imagref(op_loc), MPC_RND_IM(rnd_sin));
251 }
252
253 if (rop_cos != NULL) {
254 /* cos(-0 - i * y) = cos(+0 + i * y) = cosh(y) - i * 0,
255 cos(-0 + i * y) = cos(+0 - i * y) = cosh(y) + i * 0,
256 where y > 0 */
257 inex_cos_re = mpfr_cosh (mpc_realref (rop_cos), mpc_imagref (op_loc), MPC_RND_RE (rnd_cos));
258
259 mpfr_set_ui (mpc_imagref (rop_cos), 0ul, MPC_RND_IM (rnd_cos));
260 if (mpfr_signbit (mpc_realref (op_loc)) == mpfr_signbit (mpc_imagref (op_loc)))
261 MPFR_CHANGE_SIGN (mpc_imagref (rop_cos));
262 }
263
264 if (overlap)
265 mpc_clear (op_loc);
266
267 return MPC_INEX12 (MPC_INEX (0, inex_sin_im), MPC_INEX (inex_cos_re, 0));
268 }
269
270 /* Fix an inexact overflow, when x is +Inf or -Inf:
271 When rnd is towards zero, change x into the largest (in absolute value)
272 floating-point number.
273 Return the inexact flag. */
274 int
mpc_fix_inf(mpfr_t x,mpfr_rnd_t rnd)275 mpc_fix_inf (mpfr_t x, mpfr_rnd_t rnd)
276 {
277 MPC_ASSERT (mpfr_inf_p (x));
278 if (!MPC_IS_LIKE_RNDZ(rnd, MPFR_SIGNBIT(x)))
279 return mpfr_sgn (x);
280 else
281 {
282 if (mpfr_sgn (x) < 0)
283 mpfr_nextabove (x);
284 else
285 mpfr_nextbelow (x);
286 return -mpfr_sgn (x);
287 }
288 }
289
290 /* Fix an inexact underflow, when x is +0 or -0:
291 When rnd is away from zero, change x into the closest floating-point number.
292 Return the inexact flag. */
293 int
mpc_fix_zero(mpfr_t x,mpfr_rnd_t rnd)294 mpc_fix_zero (mpfr_t x, mpfr_rnd_t rnd)
295 {
296 if (!MPC_IS_LIKE_RNDA(rnd, MPFR_SIGNBIT(x)))
297 return mpfr_signbit (x) == 0 ? -1 : 1;
298 else
299 {
300 if (mpfr_signbit (x) == 0)
301 {
302 mpfr_nextabove (x);
303 return 1;
304 }
305 else
306 {
307 mpfr_nextbelow (x);
308 return -1;
309 }
310 }
311 }
312
313 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)314 mpc_sin_cos (mpc_ptr rop_sin, mpc_ptr rop_cos, mpc_srcptr op,
315 mpc_rnd_t rnd_sin, mpc_rnd_t rnd_cos)
316 /* Feature not documented in the texinfo file: One of rop_sin or
317 rop_cos may be NULL, in which case it is not computed, and the
318 corresponding ternary inexact value is set to 0 (exact). */
319 {
320 if (!mpc_fin_p (op))
321 return mpc_sin_cos_nonfinite (rop_sin, rop_cos, op, rnd_sin, rnd_cos);
322 else if (mpfr_zero_p (mpc_imagref (op)))
323 return mpc_sin_cos_real (rop_sin, rop_cos, op, rnd_sin, rnd_cos);
324 else if (mpfr_zero_p (mpc_realref (op)))
325 return mpc_sin_cos_imag (rop_sin, rop_cos, op, rnd_sin, rnd_cos);
326 else {
327 /* let op = a + i*b, then sin(op) = sin(a)*cosh(b) + i*cos(a)*sinh(b)
328 and cos(op) = cos(a)*cosh(b) - i*sin(a)*sinh(b).
329
330 For Re(sin(op)) (and analogously, the other parts), we use the
331 following algorithm, with rounding to nearest for all operations
332 and working precision w:
333
334 (1) x = o(sin(a))
335 (2) y = o(cosh(b))
336 (3) r = o(x*y)
337 then the error on r is at most 4 ulps, since we can write
338 r = sin(a)*cosh(b)*(1+t)^3 with |t| <= 2^(-w),
339 thus for w >= 2, r = sin(a)*cosh(b)*(1+4*t) with |t| <= 2^(-w),
340 thus the relative error is bounded by 4*2^(-w) <= 4*ulp(r).
341 */
342 mpfr_t s, c, sh, ch, sch, csh;
343 mpfr_prec_t prec;
344 int ok;
345 int inex_re, inex_im, inex_sin, inex_cos, loop = 0;
346 mpfr_exp_t saved_emin, saved_emax;
347
348 saved_emin = mpfr_get_emin ();
349 saved_emax = mpfr_get_emax ();
350 mpfr_set_emin (mpfr_get_emin_min ());
351 mpfr_set_emax (mpfr_get_emax_max ());
352
353 prec = 2;
354 if (rop_sin != NULL)
355 {
356 mp_prec_t er, ei;
357 prec = MPC_MAX (prec, MPC_MAX_PREC (rop_sin));
358 /* since the Taylor expansion of sin(x) at x=0 is x - x^3/6 + O(x^5),
359 if x <= 2^(-p), then the second term/x is about 2^(-2p)/6, thus we
360 need at least 2p+3 bits of precision. This is true only when x is
361 exactly representable in the target precision. */
362 if (MPC_MAX_PREC (op) <= prec)
363 {
364 er = mpfr_get_exp (mpc_realref (op));
365 ei = mpfr_get_exp (mpc_imagref (op));
366 /* consider the maximal exponent only */
367 er = (er < ei) ? ei : er;
368 if (er < 0)
369 if (prec < 2 * (mp_prec_t) (-er) + 3)
370 prec = 2 * (mp_prec_t) (-er) + 3;
371 }
372 }
373 if (rop_cos != NULL)
374 prec = MPC_MAX (prec, MPC_MAX_PREC (rop_cos));
375
376 mpfr_init2 (s, 2);
377 mpfr_init2 (c, 2);
378 mpfr_init2 (sh, 2);
379 mpfr_init2 (ch, 2);
380 mpfr_init2 (sch, 2);
381 mpfr_init2 (csh, 2);
382
383 do {
384 loop ++;
385 prec += (loop <= 2) ? mpc_ceil_log2 (prec) + 5 : prec / 2;
386
387 mpfr_set_prec (s, prec);
388 mpfr_set_prec (c, prec);
389 mpfr_set_prec (sh, prec);
390 mpfr_set_prec (ch, prec);
391 mpfr_set_prec (sch, prec);
392 mpfr_set_prec (csh, prec);
393
394 mpfr_sin_cos (s, c, mpc_realref(op), MPFR_RNDN);
395 mpfr_sinh_cosh (sh, ch, mpc_imagref(op), MPFR_RNDN);
396
397 ok = 1;
398
399 if (rop_sin != NULL) {
400 /* real part of sine */
401 mpfr_mul (sch, s, ch, MPFR_RNDN);
402 ok = (!mpfr_number_p (sch))
403 || mpfr_can_round (sch, prec - 2, MPFR_RNDN, MPFR_RNDZ,
404 MPC_PREC_RE (rop_sin)
405 + (MPC_RND_RE (rnd_sin) == MPFR_RNDN));
406
407 if (ok) {
408 /* imaginary part of sine */
409 mpfr_mul (csh, c, sh, MPFR_RNDN);
410 ok = (!mpfr_number_p (csh))
411 || mpfr_can_round (csh, prec - 2, MPFR_RNDN, MPFR_RNDZ,
412 MPC_PREC_IM (rop_sin)
413 + (MPC_RND_IM (rnd_sin) == MPFR_RNDN));
414 }
415 }
416
417 if (rop_cos != NULL && ok) {
418 /* real part of cosine */
419 mpfr_mul (c, c, ch, MPFR_RNDN);
420 ok = (!mpfr_number_p (c))
421 || mpfr_can_round (c, prec - 2, MPFR_RNDN, MPFR_RNDZ,
422 MPC_PREC_RE (rop_cos)
423 + (MPC_RND_RE (rnd_cos) == MPFR_RNDN));
424
425 if (ok) {
426 /* imaginary part of cosine */
427 mpfr_mul (s, s, sh, MPFR_RNDN);
428 mpfr_neg (s, s, MPFR_RNDN);
429 ok = (!mpfr_number_p (s))
430 || mpfr_can_round (s, prec - 2, MPFR_RNDN, MPFR_RNDZ,
431 MPC_PREC_IM (rop_cos)
432 + (MPC_RND_IM (rnd_cos) == MPFR_RNDN));
433 }
434 }
435 } while (ok == 0);
436
437 if (rop_sin != NULL) {
438 inex_re = mpfr_set (mpc_realref (rop_sin), sch, MPC_RND_RE (rnd_sin));
439 if (mpfr_inf_p (sch))
440 inex_re = mpc_fix_inf (mpc_realref (rop_sin), MPC_RND_RE (rnd_sin));
441 inex_im = mpfr_set (mpc_imagref (rop_sin), csh, MPC_RND_IM (rnd_sin));
442 if (mpfr_inf_p (csh))
443 inex_im = mpc_fix_inf (mpc_imagref (rop_sin), MPC_RND_IM (rnd_sin));
444 inex_sin = MPC_INEX (inex_re, inex_im);
445 }
446 else
447 inex_sin = MPC_INEX (0,0); /* return exact if not computed */
448
449 if (rop_cos != NULL) {
450 inex_re = mpfr_set (mpc_realref (rop_cos), c, MPC_RND_RE (rnd_cos));
451 if (mpfr_inf_p (c))
452 inex_re = mpc_fix_inf (mpc_realref (rop_cos), MPC_RND_RE (rnd_cos));
453 inex_im = mpfr_set (mpc_imagref (rop_cos), s, MPC_RND_IM (rnd_cos));
454 if (mpfr_inf_p (s))
455 inex_im = mpc_fix_inf (mpc_imagref (rop_cos), MPC_RND_IM (rnd_cos));
456 inex_cos = MPC_INEX (inex_re, inex_im);
457 }
458 else
459 inex_cos = MPC_INEX (0,0); /* return exact if not computed */
460
461 mpfr_clear (s);
462 mpfr_clear (c);
463 mpfr_clear (sh);
464 mpfr_clear (ch);
465 mpfr_clear (sch);
466 mpfr_clear (csh);
467
468 /* restore the exponent range, and check the range of results */
469 mpfr_set_emin (saved_emin);
470 mpfr_set_emax (saved_emax);
471 if (rop_sin != NULL)
472 {
473 inex_re = mpfr_check_range (mpc_realref (rop_sin),
474 MPC_INEX_RE (inex_sin),
475 MPC_RND_RE (rnd_sin));
476 inex_im = mpfr_check_range (mpc_imagref (rop_sin),
477 MPC_INEX_IM (inex_sin),
478 MPC_RND_IM (rnd_sin));
479 inex_sin = MPC_INEX (inex_re, inex_im);
480 }
481 if (rop_cos != NULL)
482 {
483 inex_re = mpfr_check_range (mpc_realref (rop_cos),
484 MPC_INEX_RE (inex_cos),
485 MPC_RND_RE (rnd_cos));
486 inex_im = mpfr_check_range (mpc_imagref (rop_cos),
487 MPC_INEX_IM (inex_cos),
488 MPC_RND_IM (rnd_cos));
489 inex_cos = MPC_INEX (inex_re, inex_im);
490 }
491
492 return (MPC_INEX12 (inex_sin, inex_cos));
493 }
494 }
495