1 /* mpc_pow -- Raise a complex number to the power of another complex number. 2 3 Copyright (C) 2009, 2010, 2011, 2012 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 "mpc-impl.h" 23 24 /* Return non-zero iff c+i*d is an exact square (a+i*b)^2, 25 with a, b both of the form m*2^e with m, e integers. 26 If so, returns in a+i*b the corresponding square root, with a >= 0. 27 The variables a, b must not overlap with c, d. 28 29 We have c = a^2 - b^2 and d = 2*a*b. 30 31 If one of a, b is exact, then both are (see algorithms.tex). 32 33 Case 1: a <> 0 and b <> 0. 34 Let a = m*2^e and b = n*2^f with m, e, n, f integers, m and n odd 35 (we will treat apart the case a = 0 or b = 0). 36 Then 2*a*b = m*n*2^(e+f+1), thus necessarily e+f >= -1. 37 Assume e < 0, then f >= 0, then a^2 - b^2 = m^2*2^(2e) - n^2*2^(2f) cannot 38 be an integer, since n^2*2^(2f) is an integer, and m^2*2^(2e) is not. 39 Similarly when f < 0 (and thus e >= 0). 40 Thus we have e, f >= 0, and a, b are both integers. 41 Let A = 2a^2, then eliminating b between c = a^2 - b^2 and d = 2*a*b 42 gives A^2 - 2c*A - d^2 = 0, which has solutions c +/- sqrt(c^2+d^2). 43 We thus need c^2+d^2 to be a square, and c + sqrt(c^2+d^2) --- the solution 44 we are interested in --- to be two times a square. Then b = d/(2a) is 45 necessarily an integer. 46 47 Case 2: a = 0. Then d is necessarily zero, thus it suffices to check 48 whether c = -b^2, i.e., if -c is a square. 49 50 Case 3: b = 0. Then d is necessarily zero, thus it suffices to check 51 whether c = a^2, i.e., if c is a square. 52 */ 53 static int 54 mpc_perfect_square_p (mpz_t a, mpz_t b, mpz_t c, mpz_t d) 55 { 56 if (mpz_cmp_ui (d, 0) == 0) /* case a = 0 or b = 0 */ 57 { 58 /* necessarily c < 0 here, since we have already considered the case 59 where x is real non-negative and y is real */ 60 MPC_ASSERT (mpz_cmp_ui (c, 0) < 0); 61 mpz_neg (b, c); 62 if (mpz_perfect_square_p (b)) /* case 2 above */ 63 { 64 mpz_sqrt (b, b); 65 mpz_set_ui (a, 0); 66 return 1; /* c + i*d = (0 + i*b)^2 */ 67 } 68 } 69 else /* both a and b are non-zero */ 70 { 71 if (mpz_divisible_2exp_p (d, 1) == 0) 72 return 0; /* d must be even */ 73 mpz_mul (a, c, c); 74 mpz_addmul (a, d, d); /* c^2 + d^2 */ 75 if (mpz_perfect_square_p (a)) 76 { 77 mpz_sqrt (a, a); 78 mpz_add (a, c, a); /* c + sqrt(c^2+d^2) */ 79 if (mpz_divisible_2exp_p (a, 1)) 80 { 81 mpz_tdiv_q_2exp (a, a, 1); 82 if (mpz_perfect_square_p (a)) 83 { 84 mpz_sqrt (a, a); 85 mpz_tdiv_q_2exp (b, d, 1); /* d/2 */ 86 mpz_divexact (b, b, a); /* d/(2a) */ 87 return 1; 88 } 89 } 90 } 91 } 92 return 0; /* not a square */ 93 } 94 95 /* fix the sign of Re(z) or Im(z) in case it is zero, 96 and Re(x) is zero. 97 sign_eps is 0 if Re(x) = +0, 1 if Re(x) = -0 98 sign_a is the sign bit of Im(x). 99 Assume y is an integer (does nothing otherwise). 100 */ 101 static void 102 fix_sign (mpc_ptr z, int sign_eps, int sign_a, mpfr_srcptr y) 103 { 104 int ymod4 = -1; 105 mpfr_exp_t ey; 106 mpz_t my; 107 unsigned long int t; 108 109 mpz_init (my); 110 111 ey = mpfr_get_z_exp (my, y); 112 /* normalize so that my is odd */ 113 t = mpz_scan1 (my, 0); 114 ey += (mpfr_exp_t) t; 115 mpz_tdiv_q_2exp (my, my, t); 116 /* y = my*2^ey */ 117 118 /* compute y mod 4 (in case y is an integer) */ 119 if (ey >= 2) 120 ymod4 = 0; 121 else if (ey == 1) 122 ymod4 = mpz_tstbit (my, 0) * 2; /* correct if my < 0 */ 123 else if (ey == 0) 124 { 125 ymod4 = mpz_tstbit (my, 1) * 2 + mpz_tstbit (my, 0); 126 if (mpz_cmp_ui (my , 0) < 0) 127 ymod4 = 4 - ymod4; 128 } 129 else /* y is not an integer */ 130 goto end; 131 132 if (mpfr_zero_p (mpc_realref(z))) 133 { 134 /* we assume y is always integer in that case (FIXME: prove it): 135 (eps+I*a)^y = +0 + I*a^y for y = 1 mod 4 and sign_eps = 0 136 (eps+I*a)^y = -0 - I*a^y for y = 3 mod 4 and sign_eps = 0 */ 137 MPC_ASSERT (ymod4 == 1 || ymod4 == 3); 138 if ((ymod4 == 3 && sign_eps == 0) || 139 (ymod4 == 1 && sign_eps == 1)) 140 mpfr_neg (mpc_realref(z), mpc_realref(z), GMP_RNDZ); 141 } 142 else if (mpfr_zero_p (mpc_imagref(z))) 143 { 144 /* we assume y is always integer in that case (FIXME: prove it): 145 (eps+I*a)^y = a^y - 0*I for y = 0 mod 4 and sign_a = sign_eps 146 (eps+I*a)^y = -a^y +0*I for y = 2 mod 4 and sign_a = sign_eps */ 147 MPC_ASSERT (ymod4 == 0 || ymod4 == 2); 148 if ((ymod4 == 0 && sign_a == sign_eps) || 149 (ymod4 == 2 && sign_a != sign_eps)) 150 mpfr_neg (mpc_imagref(z), mpc_imagref(z), GMP_RNDZ); 151 } 152 153 end: 154 mpz_clear (my); 155 } 156 157 /* If x^y is exactly representable (with maybe a larger precision than z), 158 round it in z and return the (mpc) inexact flag in [0, 10]. 159 160 If x^y is not exactly representable, return -1. 161 162 If intermediate computations lead to numbers of more than maxprec bits, 163 then abort and return -2 (in that case, to avoid loops, mpc_pow_exact 164 should be called again with a larger value of maxprec). 165 166 Assume one of Re(x) or Im(x) is non-zero, and y is non-zero (y is real). 167 168 Warning: z and x might be the same variable, same for Re(z) or Im(z) and y. 169 170 In case -1 or -2 is returned, z is not modified. 171 */ 172 static int 173 mpc_pow_exact (mpc_ptr z, mpc_srcptr x, mpfr_srcptr y, mpc_rnd_t rnd, 174 mpfr_prec_t maxprec) 175 { 176 mpfr_exp_t ec, ed, ey; 177 mpz_t my, a, b, c, d, u; 178 unsigned long int t; 179 int ret = -2; 180 int sign_rex = mpfr_signbit (mpc_realref(x)); 181 int sign_imx = mpfr_signbit (mpc_imagref(x)); 182 int x_imag = mpfr_zero_p (mpc_realref(x)); 183 int z_is_y = 0; 184 mpfr_t copy_of_y; 185 186 if (mpc_realref (z) == y || mpc_imagref (z) == y) 187 { 188 z_is_y = 1; 189 mpfr_init2 (copy_of_y, mpfr_get_prec (y)); 190 mpfr_set (copy_of_y, y, GMP_RNDN); 191 } 192 193 mpz_init (my); 194 mpz_init (a); 195 mpz_init (b); 196 mpz_init (c); 197 mpz_init (d); 198 mpz_init (u); 199 200 ey = mpfr_get_z_exp (my, y); 201 /* normalize so that my is odd */ 202 t = mpz_scan1 (my, 0); 203 ey += (mpfr_exp_t) t; 204 mpz_tdiv_q_2exp (my, my, t); 205 /* y = my*2^ey with my odd */ 206 207 if (x_imag) 208 { 209 mpz_set_ui (c, 0); 210 ec = 0; 211 } 212 else 213 ec = mpfr_get_z_exp (c, mpc_realref(x)); 214 if (mpfr_zero_p (mpc_imagref(x))) 215 { 216 mpz_set_ui (d, 0); 217 ed = ec; 218 } 219 else 220 { 221 ed = mpfr_get_z_exp (d, mpc_imagref(x)); 222 if (x_imag) 223 ec = ed; 224 } 225 /* x = c*2^ec + I * d*2^ed */ 226 /* equalize the exponents of x */ 227 if (ec < ed) 228 { 229 mpz_mul_2exp (d, d, (unsigned long int) (ed - ec)); 230 if ((mpfr_prec_t) mpz_sizeinbase (d, 2) > maxprec) 231 goto end; 232 } 233 else if (ed < ec) 234 { 235 mpz_mul_2exp (c, c, (unsigned long int) (ec - ed)); 236 if ((mpfr_prec_t) mpz_sizeinbase (c, 2) > maxprec) 237 goto end; 238 ec = ed; 239 } 240 /* now ec=ed and x = (c + I * d) * 2^ec */ 241 242 /* divide by two if possible */ 243 if (mpz_cmp_ui (c, 0) == 0) 244 { 245 t = mpz_scan1 (d, 0); 246 mpz_tdiv_q_2exp (d, d, t); 247 ec += (mpfr_exp_t) t; 248 } 249 else if (mpz_cmp_ui (d, 0) == 0) 250 { 251 t = mpz_scan1 (c, 0); 252 mpz_tdiv_q_2exp (c, c, t); 253 ec += (mpfr_exp_t) t; 254 } 255 else /* neither c nor d is zero */ 256 { 257 unsigned long v; 258 t = mpz_scan1 (c, 0); 259 v = mpz_scan1 (d, 0); 260 if (v < t) 261 t = v; 262 mpz_tdiv_q_2exp (c, c, t); 263 mpz_tdiv_q_2exp (d, d, t); 264 ec += (mpfr_exp_t) t; 265 } 266 267 /* now either one of c, d is odd */ 268 269 while (ey < 0) 270 { 271 /* check if x is a square */ 272 if (ec & 1) 273 { 274 mpz_mul_2exp (c, c, 1); 275 mpz_mul_2exp (d, d, 1); 276 ec --; 277 } 278 /* now ec is even */ 279 if (mpc_perfect_square_p (a, b, c, d) == 0) 280 break; 281 mpz_swap (a, c); 282 mpz_swap (b, d); 283 ec /= 2; 284 ey ++; 285 } 286 287 if (ey < 0) 288 { 289 ret = -1; /* not representable */ 290 goto end; 291 } 292 293 /* Now ey >= 0, it thus suffices to check that x^my is representable. 294 If my > 0, this is always true. If my < 0, we first try to invert 295 (c+I*d)*2^ec. 296 */ 297 if (mpz_cmp_ui (my, 0) < 0) 298 { 299 /* If my < 0, 1 / (c + I*d) = (c - I*d)/(c^2 + d^2), thus a sufficient 300 condition is that c^2 + d^2 is a power of two, assuming |c| <> |d|. 301 Assume a prime p <> 2 divides c^2 + d^2, 302 then if p does not divide c or d, 1 / (c + I*d) cannot be exact. 303 If p divides both c and d, then we can write c = p*c', d = p*d', 304 and 1 / (c + I*d) = 1/p * 1/(c' + I*d'). This shows that if 1/(c+I*d) 305 is exact, then 1/(c' + I*d') is exact too, and we are back to the 306 previous case. In conclusion, a necessary and sufficient condition 307 is that c^2 + d^2 is a power of two. 308 */ 309 /* FIXME: we could first compute c^2+d^2 mod a limb for example */ 310 mpz_mul (a, c, c); 311 mpz_addmul (a, d, d); 312 t = mpz_scan1 (a, 0); 313 if (mpz_sizeinbase (a, 2) != 1 + t) /* a is not a power of two */ 314 { 315 ret = -1; /* not representable */ 316 goto end; 317 } 318 /* replace (c,d) by (c/(c^2+d^2), -d/(c^2+d^2)) */ 319 mpz_neg (d, d); 320 ec = -ec - (mpfr_exp_t) t; 321 mpz_neg (my, my); 322 } 323 324 /* now ey >= 0 and my >= 0, and we want to compute 325 [(c + I * d) * 2^ec] ^ (my * 2^ey). 326 327 We first compute [(c + I * d) * 2^ec]^my, then square ey times. */ 328 t = mpz_sizeinbase (my, 2) - 1; 329 mpz_set (a, c); 330 mpz_set (b, d); 331 ed = ec; 332 /* invariant: (a + I*b) * 2^ed = ((c + I*d) * 2^ec)^trunc(my/2^t) */ 333 while (t-- > 0) 334 { 335 unsigned long int v, w; 336 /* square a + I*b */ 337 mpz_mul (u, a, b); 338 mpz_mul (a, a, a); 339 mpz_submul (a, b, b); 340 mpz_mul_2exp (b, u, 1); 341 ed *= 2; 342 if (mpz_tstbit (my, t)) /* multiply by c + I*d */ 343 { 344 mpz_mul (u, a, c); 345 mpz_submul (u, b, d); /* ac-bd */ 346 mpz_mul (b, b, c); 347 mpz_addmul (b, a, d); /* bc+ad */ 348 mpz_swap (a, u); 349 ed += ec; 350 } 351 /* remove powers of two in (a,b) */ 352 if (mpz_cmp_ui (a, 0) == 0) 353 { 354 w = mpz_scan1 (b, 0); 355 mpz_tdiv_q_2exp (b, b, w); 356 ed += (mpfr_exp_t) w; 357 } 358 else if (mpz_cmp_ui (b, 0) == 0) 359 { 360 w = mpz_scan1 (a, 0); 361 mpz_tdiv_q_2exp (a, a, w); 362 ed += (mpfr_exp_t) w; 363 } 364 else 365 { 366 w = mpz_scan1 (a, 0); 367 v = mpz_scan1 (b, 0); 368 if (v < w) 369 w = v; 370 mpz_tdiv_q_2exp (a, a, w); 371 mpz_tdiv_q_2exp (b, b, w); 372 ed += (mpfr_exp_t) w; 373 } 374 if ( (mpfr_prec_t) mpz_sizeinbase (a, 2) > maxprec 375 || (mpfr_prec_t) mpz_sizeinbase (b, 2) > maxprec) 376 goto end; 377 } 378 /* now a+I*b = (c+I*d)^my */ 379 380 while (ey-- > 0) 381 { 382 unsigned long sa, sb; 383 384 /* square a + I*b */ 385 mpz_mul (u, a, b); 386 mpz_mul (a, a, a); 387 mpz_submul (a, b, b); 388 mpz_mul_2exp (b, u, 1); 389 ed *= 2; 390 391 /* divide by largest 2^n possible, to avoid many loops for e.g., 392 (2+2*I)^16777216 */ 393 sa = mpz_scan1 (a, 0); 394 sb = mpz_scan1 (b, 0); 395 sa = (sa <= sb) ? sa : sb; 396 mpz_tdiv_q_2exp (a, a, sa); 397 mpz_tdiv_q_2exp (b, b, sa); 398 ed += (mpfr_exp_t) sa; 399 400 if ( (mpfr_prec_t) mpz_sizeinbase (a, 2) > maxprec 401 || (mpfr_prec_t) mpz_sizeinbase (b, 2) > maxprec) 402 goto end; 403 } 404 405 ret = mpfr_set_z (mpc_realref(z), a, MPC_RND_RE(rnd)); 406 ret = MPC_INEX(ret, mpfr_set_z (mpc_imagref(z), b, MPC_RND_IM(rnd))); 407 mpfr_mul_2si (mpc_realref(z), mpc_realref(z), ed, MPC_RND_RE(rnd)); 408 mpfr_mul_2si (mpc_imagref(z), mpc_imagref(z), ed, MPC_RND_IM(rnd)); 409 410 end: 411 mpz_clear (my); 412 mpz_clear (a); 413 mpz_clear (b); 414 mpz_clear (c); 415 mpz_clear (d); 416 mpz_clear (u); 417 418 if (ret >= 0 && x_imag) 419 fix_sign (z, sign_rex, sign_imx, (z_is_y) ? copy_of_y : y); 420 421 if (z_is_y) 422 mpfr_clear (copy_of_y); 423 424 return ret; 425 } 426 427 /* Return 1 if y*2^k is an odd integer, 0 otherwise. 428 Adapted from MPFR, file pow.c. 429 430 Examples: with k=0, check if y is an odd integer, 431 with k=1, check if y is half-an-integer, 432 with k=-1, check if y/2 is an odd integer. 433 */ 434 #define MPFR_LIMB_HIGHBIT ((mp_limb_t) 1 << (BITS_PER_MP_LIMB - 1)) 435 static int 436 is_odd (mpfr_srcptr y, mpfr_exp_t k) 437 { 438 mpfr_exp_t expo; 439 mpfr_prec_t prec; 440 mp_size_t yn; 441 mp_limb_t *yp; 442 443 expo = mpfr_get_exp (y) + k; 444 if (expo <= 0) 445 return 0; /* |y| < 1 and not 0 */ 446 447 prec = mpfr_get_prec (y); 448 if ((mpfr_prec_t) expo > prec) 449 return 0; /* y is a multiple of 2^(expo-prec), thus not odd */ 450 451 /* 0 < expo <= prec: 452 y = 1xxxxxxxxxt.zzzzzzzzzzzzzzzzzz[000] 453 expo bits (prec-expo) bits 454 455 We have to check that: 456 (a) the bit 't' is set 457 (b) all the 'z' bits are zero 458 */ 459 460 prec = ((prec - 1) / BITS_PER_MP_LIMB + 1) * BITS_PER_MP_LIMB - expo; 461 /* number of z+0 bits */ 462 463 yn = prec / BITS_PER_MP_LIMB; 464 /* yn is the index of limb containing the 't' bit */ 465 466 yp = y->_mpfr_d; 467 /* if expo is a multiple of BITS_PER_MP_LIMB, t is bit 0 */ 468 if (expo % BITS_PER_MP_LIMB == 0 ? (yp[yn] & 1) == 0 469 : yp[yn] << ((expo % BITS_PER_MP_LIMB) - 1) != MPFR_LIMB_HIGHBIT) 470 return 0; 471 while (--yn >= 0) 472 if (yp[yn] != 0) 473 return 0; 474 return 1; 475 } 476 477 /* Put in z the value of x^y, rounded according to 'rnd'. 478 Return the inexact flag in [0, 10]. */ 479 int 480 mpc_pow (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd) 481 { 482 int ret = -2, loop, x_real, x_imag, y_real, z_real = 0, z_imag = 0; 483 mpc_t t, u; 484 mpfr_prec_t p, pr, pi, maxprec; 485 int saved_underflow, saved_overflow; 486 487 /* save the underflow or overflow flags from MPFR */ 488 saved_underflow = mpfr_underflow_p (); 489 saved_overflow = mpfr_overflow_p (); 490 491 x_real = mpfr_zero_p (mpc_imagref(x)); 492 y_real = mpfr_zero_p (mpc_imagref(y)); 493 494 if (y_real && mpfr_zero_p (mpc_realref(y))) /* case y zero */ 495 { 496 if (x_real && mpfr_zero_p (mpc_realref(x))) 497 { 498 /* we define 0^0 to be (1, +0) since the real part is 499 coherent with MPFR where 0^0 gives 1, and the sign of the 500 imaginary part cannot be determined */ 501 mpc_set_ui_ui (z, 1, 0, MPC_RNDNN); 502 return 0; 503 } 504 else /* x^0 = 1 +/- i*0 even for x=NaN see algorithms.tex for the 505 sign of zero */ 506 { 507 mpfr_t n; 508 int inex, cx1; 509 int sign_zi; 510 /* cx1 < 0 if |x| < 1 511 cx1 = 0 if |x| = 1 512 cx1 > 0 if |x| > 1 513 */ 514 mpfr_init (n); 515 inex = mpc_norm (n, x, GMP_RNDN); 516 cx1 = mpfr_cmp_ui (n, 1); 517 if (cx1 == 0 && inex != 0) 518 cx1 = -inex; 519 520 sign_zi = (cx1 < 0 && mpfr_signbit (mpc_imagref (y)) == 0) 521 || (cx1 == 0 522 && mpfr_signbit (mpc_imagref (x)) != mpfr_signbit (mpc_realref (y))) 523 || (cx1 > 0 && mpfr_signbit (mpc_imagref (y))); 524 525 /* warning: mpc_set_ui_ui does not set Im(z) to -0 if Im(rnd)=RNDD */ 526 ret = mpc_set_ui_ui (z, 1, 0, rnd); 527 528 if (MPC_RND_IM (rnd) == GMP_RNDD || sign_zi) 529 mpc_conj (z, z, MPC_RNDNN); 530 531 mpfr_clear (n); 532 return ret; 533 } 534 } 535 536 if (!mpc_fin_p (x) || !mpc_fin_p (y)) 537 { 538 /* special values: exp(y*log(x)) */ 539 mpc_init2 (u, 2); 540 mpc_log (u, x, MPC_RNDNN); 541 mpc_mul (u, u, y, MPC_RNDNN); 542 ret = mpc_exp (z, u, rnd); 543 mpc_clear (u); 544 goto end; 545 } 546 547 if (x_real) /* case x real */ 548 { 549 if (mpfr_zero_p (mpc_realref(x))) /* x is zero */ 550 { 551 /* special values: exp(y*log(x)) */ 552 mpc_init2 (u, 2); 553 mpc_log (u, x, MPC_RNDNN); 554 mpc_mul (u, u, y, MPC_RNDNN); 555 ret = mpc_exp (z, u, rnd); 556 mpc_clear (u); 557 goto end; 558 } 559 560 /* Special case 1^y = 1 */ 561 if (mpfr_cmp_ui (mpc_realref(x), 1) == 0) 562 { 563 int s1, s2; 564 s1 = mpfr_signbit (mpc_realref (y)); 565 s2 = mpfr_signbit (mpc_imagref (x)); 566 567 ret = mpc_set_ui (z, +1, rnd); 568 /* the sign of the zero imaginary part is known in some cases (see 569 algorithm.tex). In such cases we have 570 (x +s*0i)^(y+/-0i) = x^y + s*sign(y)*0i 571 where s = +/-1. We extend here this rule to fix the sign of the 572 zero part. 573 574 Note that the sign must also be set explicitly when rnd=RNDD 575 because mpfr_set_ui(z_i, 0, rnd) always sets z_i to +0. 576 */ 577 if (MPC_RND_IM (rnd) == GMP_RNDD || s1 != s2) 578 mpc_conj (z, z, MPC_RNDNN); 579 goto end; 580 } 581 582 /* x^y is real when: 583 (a) x is real and y is integer 584 (b) x is real non-negative and y is real */ 585 if (y_real && (mpfr_integer_p (mpc_realref(y)) || 586 mpfr_cmp_ui (mpc_realref(x), 0) >= 0)) 587 { 588 int s1, s2; 589 s1 = mpfr_signbit (mpc_realref (y)); 590 s2 = mpfr_signbit (mpc_imagref (x)); 591 592 ret = mpfr_pow (mpc_realref(z), mpc_realref(x), mpc_realref(y), MPC_RND_RE(rnd)); 593 ret = MPC_INEX(ret, mpfr_set_ui (mpc_imagref(z), 0, MPC_RND_IM(rnd))); 594 595 /* the sign of the zero imaginary part is known in some cases 596 (see algorithm.tex). In such cases we have (x +s*0i)^(y+/-0i) 597 = x^y + s*sign(y)*0i where s = +/-1. 598 We extend here this rule to fix the sign of the zero part. 599 600 Note that the sign must also be set explicitly when rnd=RNDD 601 because mpfr_set_ui(z_i, 0, rnd) always sets z_i to +0. 602 */ 603 if (MPC_RND_IM(rnd) == GMP_RNDD || s1 != s2) 604 mpfr_neg (mpc_imagref(z), mpc_imagref(z), MPC_RND_IM(rnd)); 605 goto end; 606 } 607 608 /* (-1)^(n+I*t) is real for n integer and t real */ 609 if (mpfr_cmp_si (mpc_realref(x), -1) == 0 && mpfr_integer_p (mpc_realref(y))) 610 z_real = 1; 611 612 /* for x real, x^y is imaginary when: 613 (a) x is negative and y is half-an-integer 614 (b) x = -1 and Re(y) is half-an-integer 615 */ 616 if ((mpfr_cmp_ui (mpc_realref(x), 0) < 0) && is_odd (mpc_realref(y), 1) 617 && (y_real || mpfr_cmp_si (mpc_realref(x), -1) == 0)) 618 z_imag = 1; 619 } 620 else /* x non real */ 621 /* I^(t*I) and (-I)^(t*I) are real for t real, 622 I^(n+t*I) and (-I)^(n+t*I) are real for n even and t real, and 623 I^(n+t*I) and (-I)^(n+t*I) are imaginary for n odd and t real 624 (s*I)^n is real for n even and imaginary for n odd */ 625 if ((mpc_cmp_si_si (x, 0, 1) == 0 || mpc_cmp_si_si (x, 0, -1) == 0 || 626 (mpfr_cmp_ui (mpc_realref(x), 0) == 0 && y_real)) && 627 mpfr_integer_p (mpc_realref(y))) 628 { /* x is I or -I, and Re(y) is an integer */ 629 if (is_odd (mpc_realref(y), 0)) 630 z_imag = 1; /* Re(y) odd: z is imaginary */ 631 else 632 z_real = 1; /* Re(y) even: z is real */ 633 } 634 else /* (t+/-t*I)^(2n) is imaginary for n odd and real for n even */ 635 if (mpfr_cmpabs (mpc_realref(x), mpc_imagref(x)) == 0 && y_real && 636 mpfr_integer_p (mpc_realref(y)) && is_odd (mpc_realref(y), 0) == 0) 637 { 638 if (is_odd (mpc_realref(y), -1)) /* y/2 is odd */ 639 z_imag = 1; 640 else 641 z_real = 1; 642 } 643 644 pr = mpfr_get_prec (mpc_realref(z)); 645 pi = mpfr_get_prec (mpc_imagref(z)); 646 p = (pr > pi) ? pr : pi; 647 p += 12; /* experimentally, seems to give less than 10% of failures in 648 Ziv's strategy; probably wrong now since q is not computed */ 649 if (p < 64) 650 p = 64; 651 mpc_init2 (u, p); 652 mpc_init2 (t, p); 653 pr += MPC_RND_RE(rnd) == GMP_RNDN; 654 pi += MPC_RND_IM(rnd) == GMP_RNDN; 655 maxprec = MPC_MAX_PREC (z); 656 x_imag = mpfr_zero_p (mpc_realref(x)); 657 for (loop = 0;; loop++) 658 { 659 int ret_exp; 660 mpfr_exp_t dr, di; 661 mpfr_prec_t q; 662 663 mpc_log (t, x, MPC_RNDNN); 664 mpc_mul (t, t, y, MPC_RNDNN); 665 666 /* Compute q such that |Re (y log x)|, |Im (y log x)| < 2^q. 667 We recompute it at each loop since we might get different 668 bounds if the precision is not enough. */ 669 q = mpfr_get_exp (mpc_realref(t)) > 0 ? mpfr_get_exp (mpc_realref(t)) : 0; 670 if (mpfr_get_exp (mpc_imagref(t)) > (mpfr_exp_t) q) 671 q = mpfr_get_exp (mpc_imagref(t)); 672 673 mpfr_clear_overflow (); 674 mpfr_clear_underflow (); 675 ret_exp = mpc_exp (u, t, MPC_RNDNN); 676 if (mpfr_underflow_p () || mpfr_overflow_p ()) { 677 /* under- and overflow flags are set by mpc_exp */ 678 mpc_set (z, u, MPC_RNDNN); 679 ret = ret_exp; 680 goto exact; 681 } 682 683 /* Since the error bound is global, we have to take into account the 684 exponent difference between the real and imaginary parts. We assume 685 either the real or the imaginary part of u is not zero. 686 */ 687 dr = mpfr_zero_p (mpc_realref(u)) ? mpfr_get_exp (mpc_imagref(u)) 688 : mpfr_get_exp (mpc_realref(u)); 689 di = mpfr_zero_p (mpc_imagref(u)) ? dr : mpfr_get_exp (mpc_imagref(u)); 690 if (dr > di) 691 { 692 di = dr - di; 693 dr = 0; 694 } 695 else 696 { 697 dr = di - dr; 698 di = 0; 699 } 700 /* the term -3 takes into account the factor 4 in the complex error 701 (see algorithms.tex) plus one due to the exponent difference: if 702 z = a + I*b, where the relative error on z is at most 2^(-p), and 703 EXP(a) = EXP(b) + k, the relative error on b is at most 2^(k-p) */ 704 if ((z_imag || (p > q + 3 + dr && mpfr_can_round (mpc_realref(u), p - q - 3 - dr, GMP_RNDN, GMP_RNDZ, pr))) && 705 (z_real || (p > q + 3 + di && mpfr_can_round (mpc_imagref(u), p - q - 3 - di, GMP_RNDN, GMP_RNDZ, pi)))) 706 break; 707 708 /* if Re(u) is not known to be zero, assume it is a normal number, i.e., 709 neither zero, Inf or NaN, otherwise we might enter an infinite loop */ 710 MPC_ASSERT (z_imag || mpfr_number_p (mpc_realref(u))); 711 /* idem for Im(u) */ 712 MPC_ASSERT (z_real || mpfr_number_p (mpc_imagref(u))); 713 714 if (ret == -2) /* we did not yet call mpc_pow_exact, or it aborted 715 because intermediate computations had > maxprec bits */ 716 { 717 /* check exact cases (see algorithms.tex) */ 718 if (y_real) 719 { 720 maxprec *= 2; 721 ret = mpc_pow_exact (z, x, mpc_realref(y), rnd, maxprec); 722 if (ret != -1 && ret != -2) 723 goto exact; 724 } 725 p += dr + di + 64; 726 } 727 else 728 p += p / 2; 729 mpc_set_prec (t, p); 730 mpc_set_prec (u, p); 731 } 732 733 if (z_real) 734 { 735 /* When the result is real (see algorithm.tex for details), 736 Im(x^y) = 737 + sign(imag(y))*0i, if |x| > 1 738 + sign(imag(x))*sign(real(y))*0i, if |x| = 1 739 - sign(imag(y))*0i, if |x| < 1 740 */ 741 mpfr_t n; 742 int inex, cx1; 743 int sign_zi, sign_rex, sign_imx; 744 /* cx1 < 0 if |x| < 1 745 cx1 = 0 if |x| = 1 746 cx1 > 0 if |x| > 1 747 */ 748 749 sign_rex = mpfr_signbit (mpc_realref (x)); 750 sign_imx = mpfr_signbit (mpc_imagref (x)); 751 mpfr_init (n); 752 inex = mpc_norm (n, x, GMP_RNDN); 753 cx1 = mpfr_cmp_ui (n, 1); 754 if (cx1 == 0 && inex != 0) 755 cx1 = -inex; 756 757 sign_zi = (cx1 < 0 && mpfr_signbit (mpc_imagref (y)) == 0) 758 || (cx1 == 0 && sign_imx != mpfr_signbit (mpc_realref (y))) 759 || (cx1 > 0 && mpfr_signbit (mpc_imagref (y))); 760 761 /* copy RE(y) to n since if z==y we will destroy Re(y) below */ 762 mpfr_set_prec (n, mpfr_get_prec (mpc_realref (y))); 763 mpfr_set (n, mpc_realref (y), GMP_RNDN); 764 ret = mpfr_set (mpc_realref(z), mpc_realref(u), MPC_RND_RE(rnd)); 765 if (y_real && (x_real || x_imag)) 766 { 767 /* FIXME: with y_real we assume Im(y) is really 0, which is the case 768 for example when y comes from pow_fr, but in case Im(y) is +0 or 769 -0, we might get different results */ 770 mpfr_set_ui (mpc_imagref (z), 0, MPC_RND_IM (rnd)); 771 fix_sign (z, sign_rex, sign_imx, n); 772 ret = MPC_INEX(ret, 0); /* imaginary part is exact */ 773 } 774 else 775 { 776 ret = MPC_INEX (ret, mpfr_set_ui (mpc_imagref (z), 0, MPC_RND_IM (rnd))); 777 /* warning: mpfr_set_ui does not set Im(z) to -0 if Im(rnd) = RNDD */ 778 if (MPC_RND_IM (rnd) == GMP_RNDD || sign_zi) 779 mpc_conj (z, z, MPC_RNDNN); 780 } 781 782 mpfr_clear (n); 783 } 784 else if (z_imag) 785 { 786 ret = mpfr_set (mpc_imagref(z), mpc_imagref(u), MPC_RND_IM(rnd)); 787 /* if z is imaginary and y real, then x cannot be real */ 788 if (y_real && x_imag) 789 { 790 int sign_rex = mpfr_signbit (mpc_realref (x)); 791 792 /* If z overlaps with y we set Re(z) before checking Re(y) below, 793 but in that case y=0, which was dealt with above. */ 794 mpfr_set_ui (mpc_realref (z), 0, MPC_RND_RE (rnd)); 795 /* Note: fix_sign only does something when y is an integer, 796 then necessarily y = 1 or 3 (mod 4), and in that case the 797 sign of Im(x) is irrelevant. */ 798 fix_sign (z, sign_rex, 0, mpc_realref (y)); 799 ret = MPC_INEX(0, ret); 800 } 801 else 802 ret = MPC_INEX(mpfr_set_ui (mpc_realref(z), 0, MPC_RND_RE(rnd)), ret); 803 } 804 else 805 ret = mpc_set (z, u, rnd); 806 exact: 807 mpc_clear (t); 808 mpc_clear (u); 809 810 /* restore underflow and overflow flags from MPFR */ 811 if (saved_underflow) 812 mpfr_set_underflow (); 813 if (saved_overflow) 814 mpfr_set_overflow (); 815 816 end: 817 return ret; 818 } 819