1 /* mpc_pow -- Raise a complex number to the power of another complex number.
2
3 Copyright (C) 2009, 2010, 2011, 2012, 2014, 2015, 2016, 2018, 2020, 2022 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
mpc_perfect_square_p(mpz_t a,mpz_t b,mpz_t c,mpz_t d)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
fix_sign(mpc_ptr z,int sign_eps,int sign_a,mpfr_srcptr y)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), MPFR_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), MPFR_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
mpc_pow_exact(mpc_ptr z,mpc_srcptr x,mpfr_srcptr y,mpc_rnd_t rnd,mpfr_prec_t maxprec)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 int inex_im;
186
187 if (mpc_realref (z) == y || mpc_imagref (z) == y)
188 {
189 z_is_y = 1;
190 mpfr_init2 (copy_of_y, mpfr_get_prec (y));
191 mpfr_set (copy_of_y, y, MPFR_RNDN);
192 }
193
194 mpz_init (my);
195 mpz_init (a);
196 mpz_init (b);
197 mpz_init (c);
198 mpz_init (d);
199 mpz_init (u);
200
201 ey = mpfr_get_z_exp (my, y);
202 /* normalize so that my is odd */
203 t = mpz_scan1 (my, 0);
204 ey += (mpfr_exp_t) t;
205 mpz_tdiv_q_2exp (my, my, t);
206 /* y = my*2^ey with my odd */
207
208 if (x_imag)
209 {
210 mpz_set_ui (c, 0);
211 ec = 0;
212 }
213 else
214 ec = mpfr_get_z_exp (c, mpc_realref(x));
215 if (mpfr_zero_p (mpc_imagref(x)))
216 {
217 mpz_set_ui (d, 0);
218 ed = ec;
219 }
220 else
221 {
222 ed = mpfr_get_z_exp (d, mpc_imagref(x));
223 if (x_imag)
224 ec = ed;
225 }
226 /* x = c*2^ec + I * d*2^ed */
227 /* equalize the exponents of x */
228 if (ec < ed)
229 {
230 mpz_mul_2exp (d, d, (unsigned long int) (ed - ec));
231 if ((mpfr_prec_t) mpz_sizeinbase (d, 2) > maxprec)
232 goto end;
233 }
234 else if (ed < ec)
235 {
236 mpz_mul_2exp (c, c, (unsigned long int) (ec - ed));
237 if ((mpfr_prec_t) mpz_sizeinbase (c, 2) > maxprec)
238 goto end;
239 ec = ed;
240 }
241 /* now ec=ed and x = (c + I * d) * 2^ec */
242
243 /* divide by two if possible */
244 if (mpz_cmp_ui (c, 0) == 0)
245 {
246 t = mpz_scan1 (d, 0);
247 mpz_tdiv_q_2exp (d, d, t);
248 ec += (mpfr_exp_t) t;
249 }
250 else if (mpz_cmp_ui (d, 0) == 0)
251 {
252 t = mpz_scan1 (c, 0);
253 mpz_tdiv_q_2exp (c, c, t);
254 ec += (mpfr_exp_t) t;
255 }
256 else /* neither c nor d is zero */
257 {
258 unsigned long v;
259 t = mpz_scan1 (c, 0);
260 v = mpz_scan1 (d, 0);
261 if (v < t)
262 t = v;
263 mpz_tdiv_q_2exp (c, c, t);
264 mpz_tdiv_q_2exp (d, d, t);
265 ec += (mpfr_exp_t) t;
266 }
267
268 /* now either one of c, d is odd */
269
270 while (ey < 0)
271 {
272 /* check if x is a square */
273 if (ec & 1)
274 {
275 mpz_mul_2exp (c, c, 1);
276 mpz_mul_2exp (d, d, 1);
277 ec --;
278 }
279 /* now ec is even */
280 if (mpc_perfect_square_p (a, b, c, d) == 0)
281 break;
282 mpz_swap (a, c);
283 mpz_swap (b, d);
284 ec /= 2;
285 ey ++;
286 }
287
288 if (ey < 0)
289 {
290 ret = -1; /* not representable */
291 goto end;
292 }
293
294 /* Now ey >= 0, it thus suffices to check that x^my is representable.
295 If my > 0, this is always true. If my < 0, we first try to invert
296 (c+I*d)*2^ec.
297 */
298 if (mpz_cmp_ui (my, 0) < 0)
299 {
300 /* If my < 0, 1 / (c + I*d) = (c - I*d)/(c^2 + d^2), thus a sufficient
301 condition is that c^2 + d^2 is a power of two, assuming |c| <> |d|.
302 Assume a prime p <> 2 divides c^2 + d^2,
303 then if p does not divide c or d, 1 / (c + I*d) cannot be exact.
304 If p divides both c and d, then we can write c = p*c', d = p*d',
305 and 1 / (c + I*d) = 1/p * 1/(c' + I*d'). This shows that if 1/(c+I*d)
306 is exact, then 1/(c' + I*d') is exact too, and we are back to the
307 previous case. In conclusion, a necessary and sufficient condition
308 is that c^2 + d^2 is a power of two.
309 */
310 /* FIXME: we could first compute c^2+d^2 mod a limb for example */
311 mpz_mul (a, c, c);
312 mpz_addmul (a, d, d);
313 t = mpz_scan1 (a, 0);
314 if (mpz_sizeinbase (a, 2) != 1 + t) /* a is not a power of two */
315 {
316 ret = -1; /* not representable */
317 goto end;
318 }
319 /* replace (c,d) by (c/(c^2+d^2), -d/(c^2+d^2)) */
320 mpz_neg (d, d);
321 ec = -ec - (mpfr_exp_t) t;
322 mpz_neg (my, my);
323 }
324
325 /* now ey >= 0 and my >= 0, and we want to compute
326 [(c + I * d) * 2^ec] ^ (my * 2^ey).
327
328 We first compute [(c + I * d) * 2^ec]^my, then square ey times. */
329 t = mpz_sizeinbase (my, 2) - 1;
330 mpz_set (a, c);
331 mpz_set (b, d);
332 ed = ec;
333 /* invariant: (a + I*b) * 2^ed = ((c + I*d) * 2^ec)^trunc(my/2^t) */
334 while (t-- > 0)
335 {
336 unsigned long int v, w;
337 /* square a + I*b */
338 mpz_mul (u, a, b);
339 mpz_mul (a, a, a);
340 mpz_submul (a, b, b);
341 mpz_mul_2exp (b, u, 1);
342 ed *= 2;
343 if (mpz_tstbit (my, t)) /* multiply by c + I*d */
344 {
345 mpz_mul (u, a, c);
346 mpz_submul (u, b, d); /* ac-bd */
347 mpz_mul (b, b, c);
348 mpz_addmul (b, a, d); /* bc+ad */
349 mpz_swap (a, u);
350 ed += ec;
351 }
352 /* remove powers of two in (a,b) */
353 if (mpz_cmp_ui (a, 0) == 0)
354 {
355 w = mpz_scan1 (b, 0);
356 mpz_tdiv_q_2exp (b, b, w);
357 ed += (mpfr_exp_t) w;
358 }
359 else if (mpz_cmp_ui (b, 0) == 0)
360 {
361 w = mpz_scan1 (a, 0);
362 mpz_tdiv_q_2exp (a, a, w);
363 ed += (mpfr_exp_t) w;
364 }
365 else
366 {
367 w = mpz_scan1 (a, 0);
368 v = mpz_scan1 (b, 0);
369 if (v < w)
370 w = v;
371 mpz_tdiv_q_2exp (a, a, w);
372 mpz_tdiv_q_2exp (b, b, w);
373 ed += (mpfr_exp_t) w;
374 }
375 if ( (mpfr_prec_t) mpz_sizeinbase (a, 2) > maxprec
376 || (mpfr_prec_t) mpz_sizeinbase (b, 2) > maxprec)
377 goto end;
378 }
379 /* now a+I*b = (c+I*d)^my */
380
381 while (ey-- > 0)
382 {
383 unsigned long sa, sb;
384
385 /* square a + I*b */
386 mpz_mul (u, a, b);
387 mpz_mul (a, a, a);
388 mpz_submul (a, b, b);
389 mpz_mul_2exp (b, u, 1);
390 ed *= 2;
391
392 /* divide by largest 2^n possible, to avoid many loops for e.g.,
393 (2+2*I)^16777216 */
394 sa = mpz_scan1 (a, 0);
395 sb = mpz_scan1 (b, 0);
396 sa = (sa <= sb) ? sa : sb;
397 mpz_tdiv_q_2exp (a, a, sa);
398 mpz_tdiv_q_2exp (b, b, sa);
399 ed += (mpfr_exp_t) sa;
400
401 if ( (mpfr_prec_t) mpz_sizeinbase (a, 2) > maxprec
402 || (mpfr_prec_t) mpz_sizeinbase (b, 2) > maxprec)
403 goto end;
404 }
405
406 ret = mpfr_set_z (mpc_realref(z), a, MPC_RND_RE(rnd));
407 inex_im = mpfr_set_z (mpc_imagref(z), b, MPC_RND_IM(rnd));
408 ret = MPC_INEX(ret, inex_im);
409 mpfr_mul_2si (mpc_realref(z), mpc_realref(z), ed, MPC_RND_RE(rnd));
410 mpfr_mul_2si (mpc_imagref(z), mpc_imagref(z), ed, MPC_RND_IM(rnd));
411
412 end:
413 mpz_clear (my);
414 mpz_clear (a);
415 mpz_clear (b);
416 mpz_clear (c);
417 mpz_clear (d);
418 mpz_clear (u);
419
420 if (ret >= 0 && x_imag)
421 fix_sign (z, sign_rex, sign_imx, (z_is_y) ? copy_of_y : y);
422
423 if (z_is_y)
424 mpfr_clear (copy_of_y);
425
426 return ret;
427 }
428
429 /* Return 1 if y*2^k is an odd integer, 0 otherwise.
430 Adapted from MPFR, file pow.c.
431
432 Examples: with k=0, check if y is an odd integer,
433 with k=1, check if y is half-an-integer,
434 with k=-1, check if y/2 is an odd integer.
435 */
436 #define MPFR_LIMB_HIGHBIT ((mp_limb_t) 1 << (BITS_PER_MP_LIMB - 1))
437 static int
is_odd(mpfr_srcptr y,mpfr_exp_t k)438 is_odd (mpfr_srcptr y, mpfr_exp_t k)
439 {
440 mpfr_exp_t expo;
441 mpfr_prec_t prec;
442 mp_size_t yn;
443 mp_limb_t *yp;
444
445 expo = mpfr_get_exp (y) + k;
446 if (expo <= 0)
447 return 0; /* |y| < 1 and not 0 */
448
449 prec = mpfr_get_prec (y);
450 if ((mpfr_prec_t) expo > prec)
451 return 0; /* y is a multiple of 2^(expo-prec), thus not odd */
452
453 /* 0 < expo <= prec:
454 y = 1xxxxxxxxxt.zzzzzzzzzzzzzzzzzz[000]
455 expo bits (prec-expo) bits
456
457 We have to check that:
458 (a) the bit 't' is set
459 (b) all the 'z' bits are zero
460 */
461
462 prec = ((prec - 1) / BITS_PER_MP_LIMB + 1) * BITS_PER_MP_LIMB - expo;
463 /* number of z+0 bits */
464
465 yn = prec / BITS_PER_MP_LIMB;
466 /* yn is the index of limb containing the 't' bit */
467
468 yp = y->_mpfr_d;
469 /* if expo is a multiple of BITS_PER_MP_LIMB, t is bit 0 */
470 if (expo % BITS_PER_MP_LIMB == 0 ? (yp[yn] & 1) == 0
471 : yp[yn] << ((expo % BITS_PER_MP_LIMB) - 1) != MPFR_LIMB_HIGHBIT)
472 return 0;
473 while (--yn >= 0)
474 if (yp[yn] != 0)
475 return 0;
476 return 1;
477 }
478
479 /* Put in z the value of x^y, rounded according to 'rnd'.
480 Return the inexact flag in [0, 10]. */
481 int
mpc_pow(mpc_ptr z,mpc_srcptr x,mpc_srcptr y,mpc_rnd_t rnd)482 mpc_pow (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd)
483 {
484 int ret = -2, loop, x_real, x_imag, y_real, z_real = 0, z_imag = 0,
485 ramified = 0;
486 mpc_t t, u;
487 mpfr_prec_t p, pr, pi, maxprec;
488 int saved_underflow, saved_overflow;
489 int inex_re, inex_im;
490 mpfr_exp_t saved_emin, saved_emax;
491
492 /* save the underflow or overflow flags from MPFR */
493 saved_underflow = mpfr_underflow_p ();
494 saved_overflow = mpfr_overflow_p ();
495
496 x_real = mpfr_zero_p (mpc_imagref(x));
497 y_real = mpfr_zero_p (mpc_imagref(y));
498
499 if (y_real && mpfr_zero_p (mpc_realref(y))) /* case y zero */
500 {
501 if (x_real && mpfr_zero_p (mpc_realref(x)))
502 {
503 /* we define 0^0 to be (1, +0) since the real part is
504 coherent with MPFR where 0^0 gives 1, and the sign of the
505 imaginary part cannot be determined */
506 mpc_set_ui_ui (z, 1, 0, MPC_RNDNN);
507 return 0;
508 }
509 else /* x^0 = 1 +/- i*0 even for x=NaN see algorithms.tex for the
510 sign of zero */
511 {
512 mpfr_t n;
513 int inex, cx1;
514 int sign_zi;
515 /* cx1 < 0 if |x| < 1
516 cx1 = 0 if |x| = 1
517 cx1 > 0 if |x| > 1
518 */
519 mpfr_init (n);
520 inex = mpc_norm (n, x, MPFR_RNDN);
521 cx1 = mpfr_cmp_ui (n, 1);
522 if (cx1 == 0 && inex != 0)
523 cx1 = -inex;
524
525 sign_zi = (cx1 < 0 && mpfr_signbit (mpc_imagref (y)) == 0)
526 || (cx1 == 0
527 && mpfr_signbit (mpc_imagref (x)) != mpfr_signbit (mpc_realref (y)))
528 || (cx1 > 0 && mpfr_signbit (mpc_imagref (y)));
529
530 /* warning: mpc_set_ui_ui does not set Im(z) to -0 if Im(rnd)=RNDD */
531 ret = mpc_set_ui_ui (z, 1, 0, rnd);
532
533 if (MPC_RND_IM (rnd) == MPFR_RNDD || sign_zi)
534 mpc_conj (z, z, MPC_RNDNN);
535
536 mpfr_clear (n);
537 return ret;
538 }
539 }
540
541 if (!mpc_fin_p (x) || !mpc_fin_p (y))
542 {
543 /* special values: exp(y*log(x)) */
544 mpc_init2 (u, 2);
545 mpc_log (u, x, MPC_RNDNN);
546 mpc_mul (u, u, y, MPC_RNDNN);
547 ret = mpc_exp (z, u, rnd);
548 mpc_clear (u);
549 goto end;
550 }
551
552 if (x_real) /* case x real */
553 {
554 if (mpfr_zero_p (mpc_realref(x))) /* x is zero */
555 {
556 /* special values: exp(y*log(x)) */
557 mpc_init2 (u, 2);
558 mpc_log (u, x, MPC_RNDNN);
559 mpc_mul (u, u, y, MPC_RNDNN);
560 ret = mpc_exp (z, u, rnd);
561 mpc_clear (u);
562 goto end;
563 }
564
565 /* Special case 1^y = 1 */
566 if (mpfr_cmp_ui (mpc_realref(x), 1) == 0)
567 {
568 int s1, s2;
569 s1 = mpfr_signbit (mpc_realref (y));
570 s2 = mpfr_signbit (mpc_imagref (x));
571
572 ret = mpc_set_ui (z, +1, rnd);
573 /* the sign of the zero imaginary part is known in some cases (see
574 algorithm.tex). In such cases we have
575 (x +s*0i)^(y+/-0i) = x^y + s*sign(y)*0i
576 where s = +/-1. We extend here this rule to fix the sign of the
577 zero part.
578
579 Note that the sign must also be set explicitly when rnd=RNDD
580 because mpfr_set_ui(z_i, 0, rnd) always sets z_i to +0.
581 */
582 if (MPC_RND_IM (rnd) == MPFR_RNDD || s1 != s2)
583 mpc_conj (z, z, MPC_RNDNN);
584 goto end;
585 }
586
587 /* x^y is real when:
588 (a) x is real and y is integer
589 (b) x is real non-negative and y is real */
590 if (y_real && (mpfr_integer_p (mpc_realref(y)) ||
591 mpfr_cmp_ui (mpc_realref(x), 0) >= 0))
592 {
593 int s1, s2;
594 s1 = mpfr_signbit (mpc_realref (y));
595 s2 = mpfr_signbit (mpc_imagref (x));
596
597 ret = mpfr_pow (mpc_realref(z), mpc_realref(x), mpc_realref(y), MPC_RND_RE(rnd));
598 inex_im = mpfr_set_ui (mpc_imagref(z), 0, MPC_RND_IM(rnd));
599 ret = MPC_INEX(ret, inex_im);
600
601 /* the sign of the zero imaginary part is known in some cases
602 (see algorithm.tex). In such cases we have (x +s*0i)^(y+/-0i)
603 = x^y + s*sign(y)*0i where s = +/-1.
604 We extend here this rule to fix the sign of the zero part.
605
606 Note that the sign must also be set explicitly when rnd=RNDD
607 because mpfr_set_ui(z_i, 0, rnd) always sets z_i to +0.
608 */
609 if (MPC_RND_IM(rnd) == MPFR_RNDD || s1 != s2)
610 mpfr_neg (mpc_imagref(z), mpc_imagref(z), MPC_RND_IM(rnd));
611 goto end;
612 }
613
614 /* (-1)^(n+I*t) is real for n integer and t real */
615 if (mpfr_cmp_si (mpc_realref(x), -1) == 0 && mpfr_integer_p (mpc_realref(y)))
616 z_real = 1;
617
618 /* for x real, x^y is imaginary when:
619 (a) x is negative and y is half-an-integer
620 (b) x = -1 and Re(y) is half-an-integer
621 */
622 if ((mpfr_cmp_ui (mpc_realref(x), 0) < 0) && is_odd (mpc_realref(y), 1)
623 && (y_real || mpfr_cmp_si (mpc_realref(x), -1) == 0))
624 z_imag = 1;
625 }
626 else /* x non real */
627 /* I^(t*I) and (-I)^(t*I) are real for t real,
628 I^(n+t*I) and (-I)^(n+t*I) are real for n even and t real, and
629 I^(n+t*I) and (-I)^(n+t*I) are imaginary for n odd and t real
630 (s*I)^n is real for n even and imaginary for n odd */
631 if ((mpc_cmp_si_si (x, 0, 1) == 0 || mpc_cmp_si_si (x, 0, -1) == 0 ||
632 (mpfr_cmp_ui (mpc_realref(x), 0) == 0 && y_real)) &&
633 mpfr_integer_p (mpc_realref(y)))
634 { /* x is I or -I, and Re(y) is an integer */
635 if (is_odd (mpc_realref(y), 0))
636 z_imag = 1; /* Re(y) odd: z is imaginary */
637 else
638 z_real = 1; /* Re(y) even: z is real */
639 }
640 else /* (t+/-t*I)^(2n) is imaginary for n odd and real for n even */
641 if (mpfr_cmpabs (mpc_realref(x), mpc_imagref(x)) == 0 && y_real &&
642 mpfr_integer_p (mpc_realref(y)) && is_odd (mpc_realref(y), 0) == 0)
643 {
644 ramified = 1;
645 if (is_odd (mpc_realref(y), -1)) /* y/2 is odd */
646 z_imag = 1;
647 else
648 z_real = 1;
649 }
650
651 saved_emin = mpfr_get_emin ();
652 saved_emax = mpfr_get_emax ();
653 mpfr_set_emin (mpfr_get_emin_min ());
654 mpfr_set_emax (mpfr_get_emax_max ());
655
656 pr = mpfr_get_prec (mpc_realref(z));
657 pi = mpfr_get_prec (mpc_imagref(z));
658 p = (pr > pi) ? pr : pi;
659 p += 12; /* experimentally, seems to give less than 10% of failures in
660 Ziv's strategy; probably wrong now since q is not computed */
661 if (p < 64)
662 p = 64;
663 mpc_init2 (u, p);
664 mpc_init2 (t, p);
665 pr += MPC_RND_RE(rnd) == MPFR_RNDN;
666 pi += MPC_RND_IM(rnd) == MPFR_RNDN;
667 maxprec = MPC_MAX_PREC (z);
668 x_imag = mpfr_zero_p (mpc_realref(x));
669 for (loop = 0;; loop++)
670 {
671 int ret_exp;
672 mpfr_exp_t dr, di;
673 mpfr_prec_t q;
674
675 mpc_log (t, x, MPC_RNDNN);
676 mpc_mul (t, t, y, MPC_RNDNN);
677
678 /* Compute q such that |Re (y log x)|, |Im (y log x)| < 2^q.
679 We recompute it at each loop since we might get different
680 bounds if the precision is not enough. */
681 q = mpfr_get_exp (mpc_realref(t)) > 0 ? mpfr_get_exp (mpc_realref(t)) : 0;
682 if (mpfr_get_exp (mpc_imagref(t)) > (mpfr_exp_t) q)
683 q = mpfr_get_exp (mpc_imagref(t));
684
685 /* the signs of the real/imaginary parts of exp(t) are determined by the
686 quadrant of exp(i*imag(t)), which depends on imag(t) mod (2pi).
687 We ensure that p >= q + 64 to get enough precision, but this might
688 be not enough in corner cases (FIXME). */
689 if (p < q + 64)
690 {
691 p = q + 64;
692 goto try_again;
693 }
694
695 mpfr_clear_overflow ();
696 mpfr_clear_underflow ();
697 ret_exp = mpc_exp (u, t, MPC_RNDNN);
698 if (mpfr_underflow_p () || mpfr_overflow_p ()) {
699 /* under- and overflow flags are set by mpc_exp */
700 mpc_set (z, u, MPC_RNDNN);
701 inex_re = MPC_INEX_RE(ret_exp);
702 inex_im = MPC_INEX_IM(ret_exp);
703 if (mpfr_inf_p (mpc_realref (z)))
704 inex_re = mpc_fix_inf (mpc_realref (z), MPC_RND_RE(rnd));
705 if (mpfr_inf_p (mpc_imagref (z)))
706 {
707 if (z_real)
708 inex_im = mpfr_set_ui (mpc_imagref (z), 0, MPC_RND_IM(rnd));
709 else
710 inex_im = mpc_fix_inf (mpc_imagref (z), MPC_RND_IM(rnd));
711 }
712 ret = MPC_INEX(inex_re,inex_im);
713 goto exact;
714 }
715
716 /* Since the error bound is global, we have to take into account the
717 exponent difference between the real and imaginary parts. We assume
718 either the real or the imaginary part of u is not zero.
719 */
720 dr = mpfr_zero_p (mpc_realref(u)) ? mpfr_get_exp (mpc_imagref(u))
721 : mpfr_get_exp (mpc_realref(u));
722 di = mpfr_zero_p (mpc_imagref(u)) ? dr : mpfr_get_exp (mpc_imagref(u));
723 if (dr > di)
724 {
725 di = dr - di;
726 dr = 0;
727 }
728 else
729 {
730 dr = di - dr;
731 di = 0;
732 }
733 /* the term -3 takes into account the factor 4 in the complex error
734 (see algorithms.tex) plus one due to the exponent difference: if
735 z = a + I*b, where the relative error on z is at most 2^(-p), and
736 EXP(a) = EXP(b) + k, the relative error on b is at most 2^(k-p) */
737 if ((z_imag || (p > q + 3 + dr && mpfr_can_round (mpc_realref(u), p - q - 3 - dr, MPFR_RNDN, MPFR_RNDZ, pr))) &&
738 (z_real || (p > q + 3 + di && mpfr_can_round (mpc_imagref(u), p - q - 3 - di, MPFR_RNDN, MPFR_RNDZ, pi))))
739 break;
740
741 /* if Re(u) is not known to be zero, assume it is a normal number, i.e.,
742 neither zero, Inf or NaN, otherwise we might enter an infinite loop */
743 MPC_ASSERT (z_imag || mpfr_number_p (mpc_realref(u)));
744 /* idem for Im(u) */
745 MPC_ASSERT (z_real || mpfr_number_p (mpc_imagref(u)));
746
747 if (ret == -2) /* we did not yet call mpc_pow_exact, or it aborted
748 because intermediate computations had > maxprec bits */
749 {
750 /* check exact cases (see algorithms.tex) */
751 if (y_real)
752 {
753 maxprec *= 2;
754 ret = mpc_pow_exact (z, x, mpc_realref(y), rnd, maxprec);
755 if (ret != -1 && ret != -2)
756 goto exact;
757 }
758 p += dr + di + 64;
759 }
760 else
761 p += p / 2;
762 try_again:
763 mpc_set_prec (t, p);
764 mpc_set_prec (u, p);
765 }
766
767 if (z_real)
768 {
769 /* When the result is real (see algorithm.tex for details) and
770 x=x1 * (1 \pm i), y a positive integer divisible by 4, then
771 Im(x^y) = 0i with a sign that cannot be determined (and is thus
772 chosen as _1). Otherwise,
773 Im(x^y) =
774 + sign(imag(y))*0i, if |x| > 1
775 + sign(imag(x))*sign(real(y))*0i, if |x| = 1
776 - sign(imag(y))*0i, if |x| < 1
777 */
778 if (ramified)
779 ret = MPC_INEX (
780 mpfr_set (mpc_realref(z), mpc_realref(u), MPC_RND_RE(rnd)),
781 mpfr_set_ui (mpc_imagref (z), 0, MPFR_RNDN));
782 else {
783 mpfr_t n;
784 int inex, cx1;
785 int sign_zi, sign_rex, sign_imx;
786 /* cx1 < 0 if |x| < 1
787 cx1 = 0 if |x| = 1
788 cx1 > 0 if |x| > 1
789 */
790
791 sign_rex = mpfr_signbit (mpc_realref (x));
792 sign_imx = mpfr_signbit (mpc_imagref (x));
793 mpfr_init (n);
794 inex = mpc_norm (n, x, MPFR_RNDN);
795 cx1 = mpfr_cmp_ui (n, 1);
796 if (cx1 == 0 && inex != 0)
797 cx1 = -inex;
798
799 sign_zi = (cx1 < 0 && mpfr_signbit (mpc_imagref (y)) == 0)
800 || (cx1 == 0 && sign_imx != mpfr_signbit (mpc_realref (y)))
801 || (cx1 > 0 && mpfr_signbit (mpc_imagref (y)));
802
803 /* copy RE(y) to n since if z==y we will destroy Re(y) below */
804 mpfr_set_prec (n, mpfr_get_prec (mpc_realref (y)));
805 mpfr_set (n, mpc_realref (y), MPFR_RNDN);
806 ret = mpfr_set (mpc_realref(z), mpc_realref(u), MPC_RND_RE(rnd));
807 if (y_real && (x_real || x_imag))
808 {
809 /* FIXME: with y_real we assume Im(y) is really 0, which is the case
810 for example when y comes from pow_fr, but in case Im(y) is +0 or
811 -0, we might get different results */
812 mpfr_set_ui (mpc_imagref (z), 0, MPC_RND_IM (rnd));
813 fix_sign (z, sign_rex, sign_imx, n);
814 ret = MPC_INEX(ret, 0); /* imaginary part is exact */
815 }
816 else
817 {
818 inex_im = mpfr_set_ui (mpc_imagref (z), 0, MPC_RND_IM (rnd));
819 ret = MPC_INEX (ret, inex_im);
820 /* warning: mpfr_set_ui does not set Im(z) to -0 if Im(rnd) = RNDD */
821 if (MPC_RND_IM (rnd) == MPFR_RNDD || sign_zi)
822 mpc_conj (z, z, MPC_RNDNN);
823 }
824
825 mpfr_clear (n);
826 }
827 }
828 else if (z_imag)
829 {
830 if (ramified)
831 ret = MPC_INEX (
832 mpfr_set_ui (mpc_realref (z), 0, MPFR_RNDN),
833 mpfr_set (mpc_imagref(z), mpc_imagref(u), MPC_RND_IM(rnd)));
834 else
835 {
836 ret = mpfr_set (mpc_imagref(z), mpc_imagref(u), MPC_RND_IM(rnd));
837 /* if z is imaginary and y real, then x cannot be real */
838 if (y_real && x_imag)
839 {
840 int sign_rex = mpfr_signbit (mpc_realref (x));
841
842 /* If z overlaps with y we set Re(z) before checking Re(y) below,
843 but in that case y=0, which was dealt with above. */
844 mpfr_set_ui (mpc_realref (z), 0, MPC_RND_RE (rnd));
845 /* Note: fix_sign only does something when y is an integer,
846 then necessarily y = 1 or 3 (mod 4), and in that case the
847 sign of Im(x) is irrelevant. */
848 fix_sign (z, sign_rex, 0, mpc_realref (y));
849 ret = MPC_INEX(0, ret);
850 }
851 else
852 {
853 inex_re = mpfr_set_ui (mpc_realref(z), 0, MPC_RND_RE(rnd));
854 ret = MPC_INEX(inex_re, ret);
855 }
856 }
857 }
858 else
859 ret = mpc_set (z, u, rnd);
860 exact:
861 mpc_clear (t);
862 mpc_clear (u);
863
864 /* restore underflow and overflow flags from MPFR */
865 if (saved_underflow)
866 mpfr_set_underflow ();
867 if (saved_overflow)
868 mpfr_set_overflow ();
869
870 /* restore the exponent range, and check the range of results */
871 mpfr_set_emin (saved_emin);
872 mpfr_set_emax (saved_emax);
873 inex_re = mpfr_check_range (mpc_realref (z), MPC_INEX_RE(ret),
874 MPC_RND_RE (rnd));
875 inex_im = mpfr_check_range (mpc_imagref (z), MPC_INEX_IM(ret),
876 MPC_RND_IM (rnd));
877 ret = MPC_INEX(inex_re, inex_im);
878
879 end:
880 return ret;
881 }
882