1 /* mpc_pow -- Raise a complex number to the power of another complex number.
3 Copyright (C) 2009, 2010, 2011, 2012 INRIA
5 This file is part of GNU MPC.
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.
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
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/ .
21 #include <stdio.h> /* for MPC_ASSERT */
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.
29 We have c = a^2 - b^2 and d = 2*a*b.
31 If one of a, b is exact, then both are (see algorithms.tex).
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.
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.
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.
54 mpc_perfect_square_p (mpz_t a, mpz_t b, mpz_t c, mpz_t d)
56 if (mpz_cmp_ui (d, 0) == 0) /* case a = 0 or b = 0 */
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);
62 if (mpz_perfect_square_p (b)) /* case 2 above */
66 return 1; /* c + i*d = (0 + i*b)^2 */
69 else /* both a and b are non-zero */
71 if (mpz_divisible_2exp_p (d, 1) == 0)
72 return 0; /* d must be even */
74 mpz_addmul (a, d, d); /* c^2 + d^2 */
75 if (mpz_perfect_square_p (a))
78 mpz_add (a, c, a); /* c + sqrt(c^2+d^2) */
79 if (mpz_divisible_2exp_p (a, 1))
81 mpz_tdiv_q_2exp (a, a, 1);
82 if (mpz_perfect_square_p (a))
85 mpz_tdiv_q_2exp (b, d, 1); /* d/2 */
86 mpz_divexact (b, b, a); /* d/(2a) */
92 return 0; /* not a square */
95 /* fix the sign of Re(z) or Im(z) in case it 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).
102 fix_sign (mpc_ptr z, int sign_eps, int sign_a, mpfr_srcptr y)
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);
118 /* compute y mod 4 (in case y is an integer) */
122 ymod4 = mpz_tstbit (my, 0) * 2; /* correct if my < 0 */
125 ymod4 = mpz_tstbit (my, 1) * 2 + mpz_tstbit (my, 0);
126 if (mpz_cmp_ui (my , 0) < 0)
129 else /* y is not an integer */
132 if (mpfr_zero_p (mpc_realref(z)))
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);
142 else if (mpfr_zero_p (mpc_imagref(z)))
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);
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].
160 If x^y is not exactly representable, return -1.
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).
166 Assume one of Re(x) or Im(x) is non-zero, and y is non-zero (y is real).
168 Warning: z and x might be the same variable, same for Re(z) or Im(z) and y.
170 In case -1 or -2 is returned, z is not modified.
173 mpc_pow_exact (mpc_ptr z, mpc_srcptr x, mpfr_srcptr y, mpc_rnd_t rnd,
176 mpfr_exp_t ec, ed, ey;
177 mpz_t my, a, b, c, d, u;
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));
186 if (mpc_realref (z) == y || mpc_imagref (z) == y)
189 mpfr_init2 (copy_of_y, mpfr_get_prec (y));
190 mpfr_set (copy_of_y, y, GMP_RNDN);
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 */
213 ec = mpfr_get_z_exp (c, mpc_realref(x));
214 if (mpfr_zero_p (mpc_imagref(x)))
221 ed = mpfr_get_z_exp (d, mpc_imagref(x));
225 /* x = c*2^ec + I * d*2^ed */
226 /* equalize the exponents of x */
229 mpz_mul_2exp (d, d, (unsigned long int) (ed - ec));
230 if ((mpfr_prec_t) mpz_sizeinbase (d, 2) > maxprec)
235 mpz_mul_2exp (c, c, (unsigned long int) (ec - ed));
236 if ((mpfr_prec_t) mpz_sizeinbase (c, 2) > maxprec)
240 /* now ec=ed and x = (c + I * d) * 2^ec */
242 /* divide by two if possible */
243 if (mpz_cmp_ui (c, 0) == 0)
245 t = mpz_scan1 (d, 0);
246 mpz_tdiv_q_2exp (d, d, t);
247 ec += (mpfr_exp_t) t;
249 else if (mpz_cmp_ui (d, 0) == 0)
251 t = mpz_scan1 (c, 0);
252 mpz_tdiv_q_2exp (c, c, t);
253 ec += (mpfr_exp_t) t;
255 else /* neither c nor d is zero */
258 t = mpz_scan1 (c, 0);
259 v = mpz_scan1 (d, 0);
262 mpz_tdiv_q_2exp (c, c, t);
263 mpz_tdiv_q_2exp (d, d, t);
264 ec += (mpfr_exp_t) t;
267 /* now either one of c, d is odd */
271 /* check if x is a square */
274 mpz_mul_2exp (c, c, 1);
275 mpz_mul_2exp (d, d, 1);
279 if (mpc_perfect_square_p (a, b, c, d) == 0)
289 ret = -1; /* not representable */
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
297 if (mpz_cmp_ui (my, 0) < 0)
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.
309 /* FIXME: we could first compute c^2+d^2 mod a limb for example */
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 */
315 ret = -1; /* not representable */
318 /* replace (c,d) by (c/(c^2+d^2), -d/(c^2+d^2)) */
320 ec = -ec - (mpfr_exp_t) t;
324 /* now ey >= 0 and my >= 0, and we want to compute
325 [(c + I * d) * 2^ec] ^ (my * 2^ey).
327 We first compute [(c + I * d) * 2^ec]^my, then square ey times. */
328 t = mpz_sizeinbase (my, 2) - 1;
332 /* invariant: (a + I*b) * 2^ed = ((c + I*d) * 2^ec)^trunc(my/2^t) */
335 unsigned long int v, w;
339 mpz_submul (a, b, b);
340 mpz_mul_2exp (b, u, 1);
342 if (mpz_tstbit (my, t)) /* multiply by c + I*d */
345 mpz_submul (u, b, d); /* ac-bd */
347 mpz_addmul (b, a, d); /* bc+ad */
351 /* remove powers of two in (a,b) */
352 if (mpz_cmp_ui (a, 0) == 0)
354 w = mpz_scan1 (b, 0);
355 mpz_tdiv_q_2exp (b, b, w);
356 ed += (mpfr_exp_t) w;
358 else if (mpz_cmp_ui (b, 0) == 0)
360 w = mpz_scan1 (a, 0);
361 mpz_tdiv_q_2exp (a, a, w);
362 ed += (mpfr_exp_t) w;
366 w = mpz_scan1 (a, 0);
367 v = mpz_scan1 (b, 0);
370 mpz_tdiv_q_2exp (a, a, w);
371 mpz_tdiv_q_2exp (b, b, w);
372 ed += (mpfr_exp_t) w;
374 if ( (mpfr_prec_t) mpz_sizeinbase (a, 2) > maxprec
375 || (mpfr_prec_t) mpz_sizeinbase (b, 2) > maxprec)
378 /* now a+I*b = (c+I*d)^my */
382 unsigned long sa, sb;
387 mpz_submul (a, b, b);
388 mpz_mul_2exp (b, u, 1);
391 /* divide by largest 2^n possible, to avoid many loops for e.g.,
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;
400 if ( (mpfr_prec_t) mpz_sizeinbase (a, 2) > maxprec
401 || (mpfr_prec_t) mpz_sizeinbase (b, 2) > maxprec)
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));
418 if (ret >= 0 && x_imag)
419 fix_sign (z, sign_rex, sign_imx, (z_is_y) ? copy_of_y : y);
422 mpfr_clear (copy_of_y);
427 /* Return 1 if y*2^k is an odd integer, 0 otherwise.
428 Adapted from MPFR, file pow.c.
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.
434 #define MPFR_LIMB_HIGHBIT ((mp_limb_t) 1 << (BITS_PER_MP_LIMB - 1))
436 is_odd (mpfr_srcptr y, mpfr_exp_t k)
443 expo = mpfr_get_exp (y) + k;
445 return 0; /* |y| < 1 and not 0 */
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 */
452 y = 1xxxxxxxxxt.zzzzzzzzzzzzzzzzzz[000]
453 expo bits (prec-expo) bits
455 We have to check that:
456 (a) the bit 't' is set
457 (b) all the 'z' bits are zero
460 prec = ((prec - 1) / BITS_PER_MP_LIMB + 1) * BITS_PER_MP_LIMB - expo;
461 /* number of z+0 bits */
463 yn = prec / BITS_PER_MP_LIMB;
464 /* yn is the index of limb containing the 't' bit */
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)
477 /* Put in z the value of x^y, rounded according to 'rnd'.
478 Return the inexact flag in [0, 10]. */
480 mpc_pow (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd)
482 int ret = -2, loop, x_real, x_imag, y_real, z_real = 0, z_imag = 0;
484 mpfr_prec_t p, pr, pi, maxprec;
485 int saved_underflow, saved_overflow;
487 /* save the underflow or overflow flags from MPFR */
488 saved_underflow = mpfr_underflow_p ();
489 saved_overflow = mpfr_overflow_p ();
491 x_real = mpfr_zero_p (mpc_imagref(x));
492 y_real = mpfr_zero_p (mpc_imagref(y));
494 if (y_real && mpfr_zero_p (mpc_realref(y))) /* case y zero */
496 if (x_real && mpfr_zero_p (mpc_realref(x)))
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);
504 else /* x^0 = 1 +/- i*0 even for x=NaN see algorithms.tex for the
510 /* cx1 < 0 if |x| < 1
515 inex = mpc_norm (n, x, GMP_RNDN);
516 cx1 = mpfr_cmp_ui (n, 1);
517 if (cx1 == 0 && inex != 0)
520 sign_zi = (cx1 < 0 && mpfr_signbit (mpc_imagref (y)) == 0)
522 && mpfr_signbit (mpc_imagref (x)) != mpfr_signbit (mpc_realref (y)))
523 || (cx1 > 0 && mpfr_signbit (mpc_imagref (y)));
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);
528 if (MPC_RND_IM (rnd) == GMP_RNDD || sign_zi)
529 mpc_conj (z, z, MPC_RNDNN);
536 if (!mpc_fin_p (x) || !mpc_fin_p (y))
538 /* special values: exp(y*log(x)) */
540 mpc_log (u, x, MPC_RNDNN);
541 mpc_mul (u, u, y, MPC_RNDNN);
542 ret = mpc_exp (z, u, rnd);
547 if (x_real) /* case x real */
549 if (mpfr_zero_p (mpc_realref(x))) /* x is zero */
551 /* special values: exp(y*log(x)) */
553 mpc_log (u, x, MPC_RNDNN);
554 mpc_mul (u, u, y, MPC_RNDNN);
555 ret = mpc_exp (z, u, rnd);
560 /* Special case 1^y = 1 */
561 if (mpfr_cmp_ui (mpc_realref(x), 1) == 0)
564 s1 = mpfr_signbit (mpc_realref (y));
565 s2 = mpfr_signbit (mpc_imagref (x));
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
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.
577 if (MPC_RND_IM (rnd) == GMP_RNDD || s1 != s2)
578 mpc_conj (z, z, MPC_RNDNN);
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))
589 s1 = mpfr_signbit (mpc_realref (y));
590 s2 = mpfr_signbit (mpc_imagref (x));
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)));
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.
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.
603 if (MPC_RND_IM(rnd) == GMP_RNDD || s1 != s2)
604 mpfr_neg (mpc_imagref(z), mpc_imagref(z), MPC_RND_IM(rnd));
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)))
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
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))
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 */
632 z_real = 1; /* Re(y) even: z is real */
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)
638 if (is_odd (mpc_realref(y), -1)) /* y/2 is odd */
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 */
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++)
662 /* to avoid warning message, real initialisation below */
664 mpc_log (t, x, MPC_RNDNN);
665 mpc_mul (t, t, y, MPC_RNDNN);
668 /* compute q such that |Re (y log x)|, |Im (y log x)| < 2^q */
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));
674 mpfr_clear_overflow ();
675 mpfr_clear_underflow ();
676 ret_exp = mpc_exp (u, t, MPC_RNDNN);
677 if (mpfr_underflow_p () || mpfr_overflow_p ()) {
678 /* under- and overflow flags are set by mpc_exp */
679 mpc_set (z, u, MPC_RNDNN);
684 /* Since the error bound is global, we have to take into account the
685 exponent difference between the real and imaginary parts. We assume
686 either the real or the imaginary part of u is not zero.
688 dr = mpfr_zero_p (mpc_realref(u)) ? mpfr_get_exp (mpc_imagref(u))
689 : mpfr_get_exp (mpc_realref(u));
690 di = mpfr_zero_p (mpc_imagref(u)) ? dr : mpfr_get_exp (mpc_imagref(u));
701 /* the term -3 takes into account the factor 4 in the complex error
702 (see algorithms.tex) plus one due to the exponent difference: if
703 z = a + I*b, where the relative error on z is at most 2^(-p), and
704 EXP(a) = EXP(b) + k, the relative error on b is at most 2^(k-p) */
705 if ((z_imag || (p > q + 3 + dr && mpfr_can_round (mpc_realref(u), p - q - 3 - dr, GMP_RNDN, GMP_RNDZ, pr))) &&
706 (z_real || (p > q + 3 + di && mpfr_can_round (mpc_imagref(u), p - q - 3 - di, GMP_RNDN, GMP_RNDZ, pi))))
709 /* if Re(u) is not known to be zero, assume it is a normal number, i.e.,
710 neither zero, Inf or NaN, otherwise we might enter an infinite loop */
711 MPC_ASSERT (z_imag || mpfr_number_p (mpc_realref(u)));
713 MPC_ASSERT (z_real || mpfr_number_p (mpc_imagref(u)));
715 if (ret == -2) /* we did not yet call mpc_pow_exact, or it aborted
716 because intermediate computations had > maxprec bits */
718 /* check exact cases (see algorithms.tex) */
722 ret = mpc_pow_exact (z, x, mpc_realref(y), rnd, maxprec);
723 if (ret != -1 && ret != -2)
736 /* When the result is real (see algorithm.tex for details),
738 + sign(imag(y))*0i, if |x| > 1
739 + sign(imag(x))*sign(real(y))*0i, if |x| = 1
740 - sign(imag(y))*0i, if |x| < 1
744 int sign_zi, sign_rex, sign_imx;
745 /* cx1 < 0 if |x| < 1
750 sign_rex = mpfr_signbit (mpc_realref (x));
751 sign_imx = mpfr_signbit (mpc_imagref (x));
753 inex = mpc_norm (n, x, GMP_RNDN);
754 cx1 = mpfr_cmp_ui (n, 1);
755 if (cx1 == 0 && inex != 0)
758 sign_zi = (cx1 < 0 && mpfr_signbit (mpc_imagref (y)) == 0)
759 || (cx1 == 0 && sign_imx != mpfr_signbit (mpc_realref (y)))
760 || (cx1 > 0 && mpfr_signbit (mpc_imagref (y)));
762 /* copy RE(y) to n since if z==y we will destroy Re(y) below */
763 mpfr_set_prec (n, mpfr_get_prec (mpc_realref (y)));
764 mpfr_set (n, mpc_realref (y), GMP_RNDN);
765 ret = mpfr_set (mpc_realref(z), mpc_realref(u), MPC_RND_RE(rnd));
766 if (y_real && (x_real || x_imag))
768 /* FIXME: with y_real we assume Im(y) is really 0, which is the case
769 for example when y comes from pow_fr, but in case Im(y) is +0 or
770 -0, we might get different results */
771 mpfr_set_ui (mpc_imagref (z), 0, MPC_RND_IM (rnd));
772 fix_sign (z, sign_rex, sign_imx, n);
773 ret = MPC_INEX(ret, 0); /* imaginary part is exact */
777 ret = MPC_INEX (ret, mpfr_set_ui (mpc_imagref (z), 0, MPC_RND_IM (rnd)));
778 /* warning: mpfr_set_ui does not set Im(z) to -0 if Im(rnd) = RNDD */
779 if (MPC_RND_IM (rnd) == GMP_RNDD || sign_zi)
780 mpc_conj (z, z, MPC_RNDNN);
787 ret = mpfr_set (mpc_imagref(z), mpc_imagref(u), MPC_RND_IM(rnd));
788 /* if z is imaginary and y real, then x cannot be real */
789 if (y_real && x_imag)
791 int sign_rex = mpfr_signbit (mpc_realref (x));
793 /* If z overlaps with y we set Re(z) before checking Re(y) below,
794 but in that case y=0, which was dealt with above. */
795 mpfr_set_ui (mpc_realref (z), 0, MPC_RND_RE (rnd));
796 /* Note: fix_sign only does something when y is an integer,
797 then necessarily y = 1 or 3 (mod 4), and in that case the
798 sign of Im(x) is irrelevant. */
799 fix_sign (z, sign_rex, 0, mpc_realref (y));
800 ret = MPC_INEX(0, ret);
803 ret = MPC_INEX(mpfr_set_ui (mpc_realref(z), 0, MPC_RND_RE(rnd)), ret);
806 ret = mpc_set (z, u, rnd);
811 /* restore underflow and overflow flags from MPFR */
813 mpfr_set_underflow ();
815 mpfr_set_overflow ();