1 /*
2 Copyright (C) 2014 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 #define D_ABS(x) ((x) < 0.0 ? (-(x)) : (x))
15
16 int
acb_hypgeom_pfq_choose_n_double(slong * nn,const double * are,const double * aim,slong p,const double * bre,const double * bim,slong q,double log2_z,slong n_skip,slong n_min,slong n_max,slong prec)17 acb_hypgeom_pfq_choose_n_double(slong * nn,
18 const double * are, const double * aim, slong p,
19 const double * bre, const double * bim, slong q,
20 double log2_z,
21 slong n_skip, slong n_min, slong n_max, slong prec)
22 {
23 double increase, term, term_max, accuracy, accuracy_best, t, u;
24 double required_decrease;
25 slong k, n, n_best;
26 int success;
27
28 if (p < q)
29 required_decrease = 0.01;
30 else if (p == q)
31 required_decrease = 0.0001;
32 else
33 required_decrease = 0.01;
34
35 term = term_max = accuracy = accuracy_best = 0.0;
36 success = 0;
37
38 for (n = n_best = n_skip; n < n_max; n++)
39 {
40 t = 1.0;
41
42 for (k = 0; k < FLINT_MAX(p, q); k++)
43 {
44 if (k < p)
45 {
46 u = (are[k]+n-1)*(are[k]+n-1) + (aim[k]*aim[k]);
47 t *= D_ABS(u);
48 }
49
50 if (k < q)
51 {
52 u = (bre[k]+n-1)*(bre[k]+n-1) + (bim[k]*bim[k]);
53 u = D_ABS(u);
54
55 if (u > 1e-100)
56 t /= u;
57 }
58 }
59
60 increase = 0.5 * log(t) * 1.4426950408889634074 + log2_z;
61
62 term += increase;
63 term_max = FLINT_MAX(term_max, term);
64 accuracy = term_max - term;
65
66 if (accuracy > accuracy_best && n >= n_min && increase < -required_decrease)
67 {
68 n_best = n;
69 accuracy_best = accuracy;
70 }
71
72 if (accuracy_best > prec + 4)
73 {
74 success = 1;
75 break;
76 }
77 }
78
79 *nn = n_best;
80 return success;
81 }
82
83 slong
acb_hypgeom_pfq_choose_n_max(acb_srcptr a,slong p,acb_srcptr b,slong q,const acb_t z,slong prec,slong n_max)84 acb_hypgeom_pfq_choose_n_max(acb_srcptr a, slong p,
85 acb_srcptr b, slong q, const acb_t z,
86 slong prec, slong n_max)
87 {
88 slong n_skip, n_min, n_terminating, nint;
89 slong k, n;
90 double log2_z;
91 double * are;
92 double * aim;
93 double * bre;
94 double * bim;
95 mag_t zmag;
96 int success;
97
98 if (acb_is_zero(z) || !acb_is_finite(z))
99 return 1;
100
101 mag_init(zmag);
102
103 are = flint_malloc(sizeof(double) * 2 * (p + q));
104 aim = are + p;
105 bre = aim + p;
106 bim = bre + q;
107
108 acb_get_mag(zmag, z);
109 log2_z = mag_get_d_log2_approx(zmag);
110
111 n_skip = 1;
112 n_min = 1;
113 n_terminating = WORD_MAX;
114
115 for (k = 0; k < p; k++)
116 {
117 are[k] = arf_get_d(arb_midref(acb_realref(a + k)), ARF_RND_DOWN);
118 aim[k] = arf_get_d(arb_midref(acb_imagref(a + k)), ARF_RND_DOWN);
119
120 /* If the series is terminating, stop at this n. */
121 if (acb_is_int(a + k) && are[k] <= 0.0)
122 {
123 n_terminating = FLINT_MIN(n_terminating, (slong) (-are[k] + 1));
124 n_terminating = FLINT_MAX(n_terminating, 1);
125 }
126 else if (are[k] <= 0.01 && D_ABS(aim[k]) < 0.01)
127 {
128 /* If very close to an integer, we don't want to deal with it
129 using doubles, so fast forward (todo: could work with the
130 log2 of the difference instead when this happens). */
131 nint = floor(are[k] + 0.5);
132 if (D_ABS(nint - are[k]) < 0.01)
133 n_skip = FLINT_MAX(n_skip, 2 - nint);
134 }
135 }
136
137 n_max = FLINT_MIN(n_max, n_terminating);
138
139 for (k = 0; k < q; k++)
140 {
141 bre[k] = arf_get_d(arb_midref(acb_realref(b + k)), ARF_RND_DOWN);
142 bim[k] = arf_get_d(arb_midref(acb_imagref(b + k)), ARF_RND_DOWN);
143
144 if (bre[k] <= 0.25)
145 {
146 n_min = FLINT_MAX(n_min, 2 - bre[k]);
147
148 /* Also avoid near-integers here (can even allow exact
149 integers when computing regularized hypergeometric functions). */
150 if (bre[k] <= 0.01 && D_ABS(bim[k]) < 0.01)
151 {
152 nint = floor(bre[k] + 0.5);
153 if (D_ABS(nint - bre[k]) < 0.01)
154 n_skip = FLINT_MAX(n_skip, 2 - nint);
155 }
156 }
157 }
158
159 success = acb_hypgeom_pfq_choose_n_double(&n, are, aim, p, bre, bim, q,
160 log2_z, n_skip, n_min, n_max, prec);
161
162 if (!success)
163 {
164 if (n_terminating <= n_max)
165 {
166 n = n_terminating;
167 }
168 else
169 {
170 n = FLINT_MAX(n_min, n);
171 n = FLINT_MIN(n_max, n);
172 }
173 }
174
175 flint_free(are);
176 mag_clear(zmag);
177
178 return n;
179 }
180
181 slong
acb_hypgeom_pfq_choose_n(acb_srcptr a,slong p,acb_srcptr b,slong q,const acb_t z,slong prec)182 acb_hypgeom_pfq_choose_n(acb_srcptr a, slong p,
183 acb_srcptr b, slong q,
184 const acb_t z, slong prec)
185 {
186 return acb_hypgeom_pfq_choose_n_max(a, p, b, q, z, prec,
187 FLINT_MIN(WORD_MAX / 2, 50 + 10.0 * prec));
188 }
189
190 slong
acb_hypgeom_pfq_series_choose_n(const acb_poly_struct * a,slong p,const acb_poly_struct * b,slong q,const acb_poly_t z,slong len,slong prec)191 acb_hypgeom_pfq_series_choose_n(const acb_poly_struct * a, slong p,
192 const acb_poly_struct * b, slong q,
193 const acb_poly_t z, slong len, slong prec)
194 {
195 slong n_skip, n_min, n_max, n_terminating, nint;
196 slong k, n;
197 double log2_z;
198 double * are;
199 double * aim;
200 double * bre;
201 double * bim;
202 acb_t t;
203 mag_t zmag;
204 int success;
205
206 if (z->length == 0) /* todo: check finite */
207 return 1;
208
209 mag_init(zmag);
210 acb_init(t);
211
212 are = flint_malloc(sizeof(double) * 2 * (p + q));
213 aim = are + p;
214 bre = aim + p;
215 bim = bre + q;
216
217 acb_get_mag(zmag, z->coeffs);
218 log2_z = mag_get_d_log2_approx(zmag);
219
220 n_skip = 1;
221 n_min = 1;
222 n_max = FLINT_MIN(WORD_MAX / 2, 50 + 10.0 * prec);
223 n_terminating = WORD_MAX;
224
225 /* e.g. for exp(x + O(x^100)), make sure that we actually
226 run the recurrence up to order 100.
227 todo: compute correctly based on degrees of inputs */
228 n_min = FLINT_MAX(n_min, len);
229 n_max = FLINT_MAX(n_max, n_min);
230
231 for (k = 0; k < p; k++)
232 {
233 acb_poly_get_coeff_acb(t, a + k, 0);
234
235 are[k] = arf_get_d(arb_midref(acb_realref(t)), ARF_RND_DOWN);
236 aim[k] = arf_get_d(arb_midref(acb_imagref(t)), ARF_RND_DOWN);
237
238 /* If the series is terminating, stop at this n. */
239 if (acb_is_int(t) && are[k] <= 0.0 && acb_poly_length(a + k) <= 1)
240 {
241 n_terminating = FLINT_MIN(n_terminating, (slong) (-are[k] + 1));
242 n_terminating = FLINT_MAX(n_terminating, 1);
243 }
244 else if (are[k] <= 0.01 && D_ABS(aim[k]) < 0.01)
245 {
246 /* If very close to an integer, we don't want to deal with it
247 using doubles, so fast forward (todo: could work with the
248 log2 of the difference instead when this happens). */
249 nint = floor(are[k] + 0.5);
250 if (D_ABS(nint - are[k]) < 0.01)
251 n_skip = FLINT_MAX(n_skip, 2 - nint);
252 }
253 }
254
255 n_max = FLINT_MIN(n_max, n_terminating);
256
257 for (k = 0; k < q; k++)
258 {
259 acb_poly_get_coeff_acb(t, b + k, 0);
260
261 bre[k] = arf_get_d(arb_midref(acb_realref(t)), ARF_RND_DOWN);
262 bim[k] = arf_get_d(arb_midref(acb_imagref(t)), ARF_RND_DOWN);
263
264 if (bre[k] <= 0.25)
265 {
266 n_min = FLINT_MAX(n_min, 2 - bre[k]);
267
268 /* Also avoid near-integers here (can even allow exact
269 integers when computing regularized hypergeometric functions). */
270 if (bre[k] <= 0.01 && D_ABS(bim[k]) < 0.01)
271 {
272 nint = floor(bre[k] + 0.5);
273 if (D_ABS(nint - bre[k]) < 0.01)
274 n_skip = FLINT_MAX(n_skip, 2 - nint);
275 }
276 }
277 }
278
279 success = acb_hypgeom_pfq_choose_n_double(&n, are, aim, p, bre, bim, q,
280 log2_z, n_skip, n_min, n_max, prec);
281
282 if (!success)
283 {
284 if (n_terminating <= n_max)
285 {
286 n = n_terminating;
287 }
288 else
289 {
290 n = FLINT_MAX(n_min, n);
291 n = FLINT_MIN(n_max, n);
292 }
293 }
294
295 flint_free(are);
296 acb_clear(t);
297 mag_clear(zmag);
298
299 return n;
300 }
301
302