1 /*
2 Copyright (C) 2015 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_hypgeom.h"
13
14 static void
bsplit(acb_t A,acb_t B,acb_t C,acb_t D,const acb_t b,const acb_t z,slong n0,slong n1,slong prec)15 bsplit(acb_t A, acb_t B, acb_t C, acb_t D,
16 const acb_t b, const acb_t z, slong n0, slong n1, slong prec)
17 {
18 if (n1 - n0 == 1)
19 {
20 acb_zero(A);
21 acb_one(B);
22 acb_neg(C, b);
23 acb_add_si(C, C, 2 - n0, prec);
24 acb_mul_si(C, C, n0 - 1, prec);
25 acb_sub(D, z, b, prec);
26 acb_add_si(D, D, 2 - 2 * n0, prec);
27 }
28 else
29 {
30 slong m;
31 acb_t T, A2, B2, C2, D2;
32
33 acb_init(T);
34 acb_init(A2);
35 acb_init(B2);
36 acb_init(C2);
37 acb_init(D2);
38
39 m = n0 + (n1 - n0) / 2;
40
41 bsplit(A, B, C, D, b, z, n0, m, prec);
42 bsplit(A2, B2, C2, D2, b, z, m, n1, prec);
43
44 acb_set(T, A);
45 acb_mul(A, A, A2, prec);
46 acb_addmul(A, B2, C, prec);
47 acb_mul(C, C, D2, prec);
48 acb_addmul(C, C2, T, prec);
49
50 acb_set(T, B);
51 acb_mul(B, B, A2, prec);
52 acb_addmul(B, B2, D, prec);
53 acb_mul(D, D, D2, prec);
54 acb_addmul(D, C2, T, prec);
55
56 acb_clear(T);
57 acb_clear(A2);
58 acb_clear(B2);
59 acb_clear(C2);
60 acb_clear(D2);
61 }
62 }
63
64 void
acb_hypgeom_u_si_rec(acb_t res,slong a,const acb_t b,const acb_t z,slong prec)65 acb_hypgeom_u_si_rec(acb_t res, slong a, const acb_t b, const acb_t z, slong prec)
66 {
67 slong k;
68 acb_t u0, u1, t;
69
70 if (a > 0)
71 flint_abort();
72
73 if (a == 0)
74 {
75 acb_one(res);
76 return;
77 }
78 else if (a == -1)
79 {
80 acb_sub(res, z, b, prec);
81 return;
82 }
83
84 /* special-case U(-n, -n+1, z) = z^n */
85 if (acb_equal_si(b, a + 1))
86 {
87 acb_pow_si(res, z, -a, prec);
88 return;
89 }
90
91 acb_init(u0);
92 acb_init(u1);
93 acb_init(t);
94
95 acb_one(u0);
96 acb_sub(u1, z, b, prec);
97
98 if (-a < 20)
99 {
100 for (k = 2; k <= -a; k++)
101 {
102 acb_neg(t, b);
103 acb_add_si(t, t, 2 - k, prec);
104 acb_mul_si(t, t, k - 1, prec);
105 acb_mul(u0, u0, t, prec);
106
107 acb_sub(t, z, b, prec);
108 acb_add_si(t, t, 2 - 2 * k, prec);
109 acb_addmul(u0, u1, t, prec);
110
111 acb_swap(u0, u1);
112 }
113
114 acb_set(res, u1);
115 }
116 else
117 {
118 acb_t A, B, C, D;
119
120 acb_init(A);
121 acb_init(B);
122 acb_init(C);
123 acb_init(D);
124
125 bsplit(A, B, C, D, b, z, 2, -a + 1, prec);
126
127 acb_sub(A, z, b, prec);
128 acb_mul(D, D, A, prec);
129 acb_add(res, C, D, prec);
130
131 acb_clear(A);
132 acb_clear(B);
133 acb_clear(C);
134 acb_clear(D);
135 }
136
137 acb_clear(u0);
138 acb_clear(u1);
139 acb_clear(t);
140 }
141
142 void
acb_hypgeom_u_1f1_series(acb_poly_t res,const acb_poly_t a,const acb_poly_t b,const acb_poly_t z,slong len,slong prec)143 acb_hypgeom_u_1f1_series(acb_poly_t res,
144 const acb_poly_t a, const acb_poly_t b, const acb_poly_t z,
145 slong len, slong prec)
146 {
147 acb_poly_t s, u, A, B;
148 acb_poly_struct aa[3];
149 arb_t c;
150 slong wlen;
151 int singular;
152
153 acb_poly_init(s);
154 acb_poly_init(u);
155 acb_poly_init(A);
156 acb_poly_init(B);
157 acb_poly_init(aa + 0);
158 acb_poly_init(aa + 1);
159 acb_poly_init(aa + 2);
160 arb_init(c);
161
162 singular = (b->length == 0) || acb_is_int(b->coeffs);
163 wlen = len + (singular != 0);
164
165 /* A = rgamma(a-b+1) * 1F~1(a,b,z) */
166 acb_poly_sub(u, a, b, prec);
167 acb_poly_add_si(u, u, 1, prec);
168 acb_poly_rgamma_series(A, u, wlen, prec);
169
170 /* todo: handle a = 1 efficiently */
171 acb_poly_set(aa, a);
172 acb_poly_set(aa + 1, b);
173 acb_poly_one(aa + 2);
174 acb_hypgeom_pfq_series_direct(s, aa, 1, aa + 1, 2, z, 1, -1, wlen, prec);
175 acb_poly_mullow(A, A, s, wlen, prec);
176
177 /* B = rgamma(a) * 1F~1(a-b+1,2-b,z) * z^(1-b) */
178 acb_poly_set(aa, u);
179 acb_poly_add_si(aa + 1, b, -2, prec);
180 acb_poly_neg(aa + 1, aa + 1);
181 acb_hypgeom_pfq_series_direct(s, aa, 1, aa + 1, 2, z, 1, -1, wlen, prec);
182 acb_poly_rgamma_series(B, a, wlen, prec);
183 acb_poly_mullow(B, B, s, wlen, prec);
184
185 acb_poly_add_si(u, b, -1, prec);
186 acb_poly_neg(u, u);
187 acb_poly_pow_series(s, z, u, wlen, prec);
188 acb_poly_mullow(B, B, s, wlen, prec);
189
190 acb_poly_sub(A, A, B, prec);
191
192 /* multiply by pi csc(pi b) */
193 acb_poly_sin_pi_series(B, b, wlen, prec);
194
195 if (singular)
196 {
197 acb_poly_shift_right(A, A, 1);
198 acb_poly_shift_right(B, B, 1);
199 }
200
201 acb_poly_div_series(res, A, B, len, prec);
202
203 arb_const_pi(c, prec);
204 _acb_vec_scalar_mul_arb(res->coeffs, res->coeffs, res->length, c, prec);
205
206 acb_poly_clear(s);
207 acb_poly_clear(u);
208 acb_poly_clear(A);
209 acb_poly_clear(B);
210 acb_poly_clear(aa + 0);
211 acb_poly_clear(aa + 1);
212 acb_poly_clear(aa + 2);
213 arb_clear(c);
214 }
215
216 void
acb_hypgeom_u_1f1(acb_t res,const acb_t a,const acb_t b,const acb_t z,slong prec)217 acb_hypgeom_u_1f1(acb_t res, const acb_t a, const acb_t b, const acb_t z, slong prec)
218 {
219 if (acb_is_int(b))
220 {
221 acb_poly_t aa, bb, zz;
222
223 acb_poly_init(aa);
224 acb_poly_init(bb);
225 acb_poly_init(zz);
226
227 acb_poly_set_acb(aa, a);
228 acb_poly_set_coeff_acb(bb, 0, b);
229 acb_poly_set_coeff_si(bb, 1, 1);
230 acb_poly_set_acb(zz, z);
231
232 acb_hypgeom_u_1f1_series(zz, aa, bb, zz, 1, prec);
233 acb_poly_get_coeff_acb(res, zz, 0);
234
235 acb_poly_clear(aa);
236 acb_poly_clear(bb);
237 acb_poly_clear(zz);
238 }
239 else
240 {
241 acb_t t, u, v;
242 acb_struct aa[3];
243
244 acb_init(t);
245 acb_init(u);
246 acb_init(v);
247 acb_init(aa + 0);
248 acb_init(aa + 1);
249 acb_init(aa + 2);
250
251 acb_set(aa, a);
252 acb_set(aa + 1, b);
253 acb_one(aa + 2);
254 acb_hypgeom_pfq_direct(u, aa, 1, aa + 1, 2, z, -1, prec);
255 acb_sub(aa, a, b, prec);
256 acb_add_ui(aa, aa, 1, prec);
257 acb_sub_ui(aa + 1, b, 2, prec);
258 acb_neg(aa + 1, aa + 1);
259 acb_hypgeom_pfq_direct(v, aa, 1, aa + 1, 2, z, -1, prec);
260
261 acb_sub_ui(aa + 1, b, 1, prec);
262
263 /* rgamma(a-b+1) * gamma(1-b) * u */
264 acb_rgamma(t, aa, prec);
265 acb_mul(u, u, t, prec);
266 acb_neg(t, aa + 1);
267 acb_gamma(t, t, prec);
268 acb_mul(u, u, t, prec);
269
270 /* rgamma(a) * gamma(b-1) * z^(1-b) * v */
271 acb_rgamma(t, a, prec);
272 acb_mul(v, v, t, prec);
273 acb_gamma(t, aa + 1, prec);
274 acb_mul(v, v, t, prec);
275 acb_neg(t, aa + 1);
276 acb_pow(t, z, t, prec);
277 acb_mul(v, v, t, prec);
278
279 acb_add(res, u, v, prec);
280
281 acb_clear(t);
282 acb_clear(u);
283 acb_clear(v);
284 acb_clear(aa + 0);
285 acb_clear(aa + 1);
286 acb_clear(aa + 2);
287 }
288 }
289
290 void
acb_hypgeom_u_choose(int * asymp,slong * wp,const acb_t a,const acb_t b,const acb_t z,slong prec)291 acb_hypgeom_u_choose(int * asymp, slong * wp,
292 const acb_t a, const acb_t b, const acb_t z, slong prec)
293 {
294 double x, y, t, cancellation;
295 double input_accuracy, direct_accuracy, asymp_accuracy;
296
297 *asymp = 0;
298 *wp = prec;
299
300 input_accuracy = acb_rel_one_accuracy_bits(z);
301 t = acb_rel_one_accuracy_bits(a);
302 input_accuracy = FLINT_MIN(input_accuracy, t);
303 t = acb_rel_one_accuracy_bits(b);
304 input_accuracy = FLINT_MIN(input_accuracy, t);
305 input_accuracy = FLINT_MAX(input_accuracy, 0.0);
306
307 /* From here we ignore the values of a, b. Taking them into account is
308 a possible future improvement... */
309
310 /* Tiny |z|. */
311 if ((arf_cmpabs_2exp_si(arb_midref(acb_realref(z)), 2) < 0 &&
312 arf_cmpabs_2exp_si(arb_midref(acb_imagref(z)), 2) < 0))
313 {
314 *asymp = 0;
315 *wp = FLINT_MAX(2, FLINT_MIN(input_accuracy + 20, prec));
316 return;
317 }
318
319 /* Huge |z|. */
320 if ((arf_cmpabs_2exp_si(arb_midref(acb_realref(z)), 64) > 0 ||
321 arf_cmpabs_2exp_si(arb_midref(acb_imagref(z)), 64) > 0))
322 {
323 *asymp = 1;
324 *wp = FLINT_MAX(2, FLINT_MIN(input_accuracy + 20, prec));
325 return;
326 }
327
328 x = arf_get_d(arb_midref(acb_realref(z)), ARF_RND_DOWN);
329 y = arf_get_d(arb_midref(acb_imagref(z)), ARF_RND_DOWN);
330
331 asymp_accuracy = sqrt(x * x + y * y) * 1.44269504088896;
332
333 /* The Kummer transformation gives less cancellation with the 1F1 series.
334 if (x < 0.0)
335 {
336 *kummer = 1;
337 x = -x;
338 } */
339
340 if (asymp_accuracy >= prec)
341 {
342 *asymp = 1;
343 *wp = FLINT_MAX(2, FLINT_MIN(input_accuracy + 20, prec));
344 return;
345 }
346
347 /* Assume U ~ 1, M ~ exp(|z|) (there is cancellation both in the
348 evaluation of M and in the linear combination) -- a better estimate
349 would account for a, b. */
350 cancellation = sqrt(x * x + y * y) * 1.44269504088896 + 5;
351
352 direct_accuracy = input_accuracy - cancellation;
353
354 if (direct_accuracy > asymp_accuracy)
355 {
356 *asymp = 0;
357 *wp = FLINT_MAX(2, FLINT_MIN(input_accuracy + 20, prec + cancellation));
358 }
359 else
360 {
361 *asymp = 1;
362 *wp = FLINT_MAX(2, FLINT_MIN(input_accuracy + 20, prec));
363 }
364 }
365
366 void
acb_hypgeom_u(acb_t res,const acb_t a,const acb_t b,const acb_t z,slong prec)367 acb_hypgeom_u(acb_t res, const acb_t a, const acb_t b, const acb_t z, slong prec)
368 {
369 acb_t t;
370 arf_srcptr av, tv;
371
372 av = arb_midref(acb_realref(a));
373
374 /* Handle small polynomial cases without divisions. */
375 /* todo: should incorporate a -> 1+a-b transformation, also... */
376 if (acb_is_int(a) && arf_sgn(av) <= 0)
377 {
378 if (arf_cmpabs_ui(av, 30) < 0 ||
379 (arf_cmpabs_ui(av, prec) < 0 &&
380 ((acb_bits(b) < 0.1 * prec && acb_bits(z) < 0.1 * prec)
381 || acb_contains_zero(z))))
382 {
383 acb_hypgeom_u_si_rec(res, arf_get_si(av, ARF_RND_DOWN), b, z, prec);
384 return;
385 }
386 }
387
388 acb_init(t);
389 acb_sub(t, a, b, prec);
390 acb_add_ui(t, t, 1, prec);
391 tv = arb_midref(acb_realref(t));
392
393 /* todo: combine these conditions with the code below */
394 if ((acb_is_int(a) && arf_sgn(av) <= 0) ||
395 (acb_is_int(t) && arf_sgn(tv) <= 0) ||
396 acb_hypgeom_u_use_asymp(z, prec))
397 {
398 acb_neg(t, a);
399 acb_pow(t, z, t, prec);
400 acb_hypgeom_u_asymp(res, a, b, z, -1, prec);
401 acb_mul(res, res, t, prec);
402 }
403 else
404 {
405 slong wp;
406 int asymp;
407
408 acb_hypgeom_u_choose(&asymp, &wp, a, b, z, prec);
409
410 if (asymp)
411 {
412 acb_neg(t, a);
413 acb_pow(t, z, t, prec);
414 acb_hypgeom_u_asymp(res, a, b, z, -1, wp);
415 acb_mul(res, res, t, prec);
416 }
417 else
418 {
419 acb_hypgeom_u_1f1(res, a, b, z, wp);
420 acb_set_round(res, res, prec);
421 }
422 }
423
424 acb_clear(t);
425 }
426
427