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 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 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 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 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