1 /*
2 Copyright (C) 2012 Fredrik Johansson
3
4 This file is part of Arb.
5
6 Arb is free software: you can redistribute it and/or modify it under
7 the terms of the GNU Lesser General Public License (LGPL) as published
8 by the Free Software Foundation; either version 2.1 of the License, or
9 (at your option) any later version. See <http://www.gnu.org/licenses/>.
10 */
11
12 #include "acb.h"
13
14 void
acb_pow_fmpz_binexp(acb_t y,const acb_t b,const fmpz_t e,slong prec)15 acb_pow_fmpz_binexp(acb_t y, const acb_t b, const fmpz_t e, slong prec)
16 {
17 slong i, wp, bits;
18
19 if (-WORD(2) <= *e && *e <= WORD(4))
20 {
21 if (*e == WORD(0))
22 {
23 acb_one(y);
24 }
25 else if (*e == WORD(1))
26 {
27 acb_set_round(y, b, prec);
28 }
29 else if (*e == -WORD(1))
30 {
31 acb_inv(y, b, prec);
32 }
33 else if (*e == WORD(2))
34 {
35 acb_mul(y, b, b, prec);
36 }
37 else if (*e == WORD(3))
38 {
39 acb_cube(y, b, prec);
40 }
41 else if (*e == WORD(4))
42 {
43 acb_mul(y, b, b, prec);
44 acb_mul(y, y, y, prec);
45 }
46 else
47 {
48 acb_inv(y, b, prec);
49 acb_mul(y, y, y, prec);
50 }
51 return;
52 }
53
54 if (fmpz_sgn(e) < 0)
55 {
56 fmpz_t f;
57 fmpz_init(f);
58 fmpz_neg(f, e);
59
60 if (acb_is_exact(b))
61 {
62 acb_pow_fmpz_binexp(y, b, f, prec + 2);
63 acb_inv(y, y, prec);
64 }
65 else
66 {
67 acb_inv(y, b, prec + fmpz_bits(e) + 2);
68 acb_pow_fmpz_binexp(y, y, f, prec);
69 }
70
71 fmpz_clear(f);
72 return;
73 }
74
75 if (!COEFF_IS_MPZ(*e) && ((*e) % 3 == 0))
76 {
77 fmpz e3 = (*e) / 3;
78 acb_pow_fmpz_binexp(y, b, &e3, prec);
79 acb_cube(y, y, prec);
80 return;
81 }
82
83 if (y == b)
84 {
85 acb_t t;
86 acb_init(t);
87 acb_set(t, b);
88 acb_pow_fmpz_binexp(y, t, e, prec);
89 acb_clear(t);
90 return;
91 }
92
93 acb_set(y, b);
94
95 bits = fmpz_bits(e);
96 wp = ARF_PREC_ADD(prec, bits);
97
98 for (i = bits - 2; i >= 0; i--)
99 {
100 acb_mul(y, y, y, wp);
101 if (fmpz_tstbit(e, i))
102 acb_mul(y, y, b, wp);
103 }
104 }
105
106 void
acb_pow_fmpz(acb_t y,const acb_t b,const fmpz_t e,slong prec)107 acb_pow_fmpz(acb_t y, const acb_t b, const fmpz_t e, slong prec)
108 {
109 acb_pow_fmpz_binexp(y, b, e, prec);
110 }
111
112 void
acb_pow_ui(acb_t y,const acb_t b,ulong e,slong prec)113 acb_pow_ui(acb_t y, const acb_t b, ulong e, slong prec)
114 {
115 fmpz_t f;
116 fmpz_init_set_ui(f, e);
117 acb_pow_fmpz(y, b, f, prec);
118 fmpz_clear(f);
119 }
120
121 void
acb_pow_si(acb_t y,const acb_t b,slong e,slong prec)122 acb_pow_si(acb_t y, const acb_t b, slong e, slong prec)
123 {
124 fmpz_t f;
125 fmpz_init(f);
126 fmpz_set_si(f, e);
127 acb_pow_fmpz(y, b, f, prec);
128 fmpz_clear(f);
129 }
130
131 void
_acb_pow_exp(acb_t z,const acb_t x,const acb_t y,slong prec)132 _acb_pow_exp(acb_t z, const acb_t x, const acb_t y, slong prec)
133 {
134 acb_t t;
135 acb_init(t);
136 acb_log(t, x, prec);
137 acb_mul(t, t, y, prec);
138 acb_exp(z, t, prec);
139 acb_clear(t);
140 }
141
142 void
_acb_pow_arb_exp(acb_t z,const acb_t x,const arb_t y,slong prec)143 _acb_pow_arb_exp(acb_t z, const acb_t x, const arb_t y, slong prec)
144 {
145 acb_t t;
146 acb_init(t);
147 acb_log(t, x, prec);
148 acb_mul_arb(t, t, y, prec);
149 acb_exp(z, t, prec);
150 acb_clear(t);
151 }
152
153 #define BINEXP_LIMIT 64
154
155 void
acb_pow_arb(acb_t z,const acb_t x,const arb_t y,slong prec)156 acb_pow_arb(acb_t z, const acb_t x, const arb_t y, slong prec)
157 {
158 const arf_struct * ymid = arb_midref(y);
159 const mag_struct * yrad = arb_radref(y);
160
161 if (arb_is_zero(y))
162 {
163 acb_one(z);
164 return;
165 }
166
167 if (acb_is_zero(x))
168 {
169 if (arb_is_positive(y))
170 acb_zero(z);
171 else
172 acb_indeterminate(z);
173 return;
174 }
175
176 if (mag_is_zero(yrad))
177 {
178 /* small half-integer or integer */
179 if (arf_cmpabs_2exp_si(ymid, BINEXP_LIMIT) < 0 &&
180 arf_is_int_2exp_si(ymid, -1))
181 {
182 fmpz_t e;
183 fmpz_init(e);
184
185 if (arf_is_int(ymid))
186 {
187 arf_get_fmpz_fixed_si(e, ymid, 0);
188 acb_pow_fmpz_binexp(z, x, e, prec);
189 }
190 else
191 {
192 arf_get_fmpz_fixed_si(e, ymid, -1);
193 if (fmpz_sgn(e) >= 0)
194 {
195 acb_sqrt(z, x, prec + fmpz_bits(e));
196 acb_pow_fmpz_binexp(z, z, e, prec);
197 }
198 else
199 {
200 fmpz_neg(e, e);
201 acb_rsqrt(z, x, prec + fmpz_bits(e));
202 acb_pow_fmpz_binexp(z, z, e, prec);
203 }
204 }
205
206 fmpz_clear(e);
207 return;
208 }
209 }
210
211 _acb_pow_arb_exp(z, x, y, prec);
212 }
213
214 void
acb_pow(acb_t z,const acb_t x,const acb_t y,slong prec)215 acb_pow(acb_t z, const acb_t x, const acb_t y, slong prec)
216 {
217 if (arb_is_zero(acb_imagref(y)))
218 {
219 acb_pow_arb(z, x, acb_realref(y), prec);
220 }
221 else
222 {
223 if (acb_is_zero(x))
224 {
225 if (arb_is_positive(acb_realref(y)))
226 acb_zero(z);
227 else
228 acb_indeterminate(z);
229 return;
230 }
231
232 _acb_pow_exp(z, x, y, prec);
233 }
234 }
235
236 void
acb_pow_analytic(acb_ptr res,const acb_t z,const acb_t w,int analytic,slong prec)237 acb_pow_analytic(acb_ptr res, const acb_t z, const acb_t w, int analytic, slong prec)
238 {
239 if (analytic && !acb_is_int(w) && arb_contains_zero(acb_imagref(z)) &&
240 !arb_is_positive(acb_realref(z)))
241 {
242 acb_indeterminate(res);
243 }
244 else
245 {
246 acb_pow(res, z, w, prec);
247 }
248 }
249
250