1185029Spjd/* mpc_sin_cos -- combined sine and cosine of a complex number. 2185029Spjd 3185029SpjdCopyright (C) 2010, 2011, 2012, 2020 INRIA 4185029Spjd 5185029SpjdThis file is part of GNU MPC. 6185029Spjd 7185029SpjdGNU MPC is free software; you can redistribute it and/or modify it under 8185029Spjdthe terms of the GNU Lesser General Public License as published by the 9185029SpjdFree Software Foundation; either version 3 of the License, or (at your 10185029Spjdoption) any later version. 11185029Spjd 12185029SpjdGNU MPC is distributed in the hope that it will be useful, but WITHOUT ANY 13185029SpjdWARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 14185029SpjdFOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for 15185029Spjdmore details. 16185029Spjd 17185029SpjdYou should have received a copy of the GNU Lesser General Public License 18185029Spjdalong with this program. If not, see http://www.gnu.org/licenses/ . 19185029Spjd*/ 20185029Spjd 21185029Spjd#include <stdio.h> 22185029Spjd#include "mpc-impl.h" 23185029Spjd 24185029Spjdstatic int 25185029Spjdmpc_sin_cos_nonfinite (mpc_ptr rop_sin, mpc_ptr rop_cos, mpc_srcptr op, 26185029Spjd mpc_rnd_t rnd_sin, mpc_rnd_t rnd_cos) 27185029Spjd /* assumes that op (that is, its real or imaginary part) is not finite */ 28185029Spjd{ 29185029Spjd int overlap; 30185029Spjd 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 178static int 179mpc_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 229static int 230mpc_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. */ 274int 275mpc_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. */ 293int 294mpc_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 313int 314mpc_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