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