1 /*
2 Copyright (C) 2017 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 "arb_hypgeom.h"
13
14 /* todo: improve for small k */
log2_bin_uiui_fast(ulong n,ulong k)15 static double log2_bin_uiui_fast(ulong n, ulong k)
16 {
17 static const float htab[] = {0.2007, 0.3374, 0.4490, 0.5437, 0.6254, 0.6963,
18 0.7580, 0.8114, 0.8572, 0.8961, 0.9285, 0.9545, 0.9746,
19 0.9888, 0.9973, 1.0};
20
21 if (k == 0 || k >= n)
22 return 0;
23 if (k > n / 2)
24 k = n - k;
25 k = (32.0 * k) / n;
26 return n * htab[FLINT_MIN(k, 15)];
27 }
28
29 /* x is unused but part of the API */
30 void
arb_hypgeom_legendre_p_ui_deriv_bound(mag_t dp,mag_t dp2,ulong n,const arb_t x,const arb_t x2sub1)31 arb_hypgeom_legendre_p_ui_deriv_bound(mag_t dp, mag_t dp2, ulong n, const arb_t x, const arb_t x2sub1)
32 {
33 mag_t c, t, u;
34
35 mag_init(c);
36 mag_init(t);
37 mag_init(u);
38
39 arb_get_mag_lower(t, x2sub1);
40 mag_rsqrt(t, t); /* t >= 1/(1-x^2)^(1/2) */
41 mag_mul_ui(u, t, n); /* u >= n/(1-x^2)^(1/2) */
42 mag_set_ui_2exp_si(c, 409, -8); /* c >= 2^(3/2)/sqrt(pi) */
43 mag_sqrt(dp, u);
44 mag_mul(dp, dp, t);
45 mag_mul(dp, dp, c); /* dp >= c*sqrt(n)/(1-x^2)^(3/4) */
46
47 mag_mul(dp2, dp, u);
48 mag_mul_2exp_si(dp2, dp2, 1); /* dp2 >= 2*c*n^(3/2)/(1-x^2)^(5/4) */
49
50 mag_set_ui(t, n);
51 mag_add_ui(t, t, 1);
52 mag_mul(t, t, t); /* t >= (n+1)^2 */
53 mag_mul_2exp_si(u, t, -1); /* u >= (n+1)^2/2 */
54 mag_min(dp, dp, u); /* |P'(x)| <= dp */
55
56 mag_mul(t, t, t);
57 mag_mul_2exp_si(u, t, -3);
58 mag_min(dp2, dp2, u); /* |P''(x)| <= dp2 */
59
60 mag_clear(c);
61 mag_clear(t);
62 mag_clear(u);
63 }
64
65 void
arb_hypgeom_legendre_p_ui(arb_t res,arb_t res_prime,ulong n,const arb_t x,slong prec)66 arb_hypgeom_legendre_p_ui(arb_t res, arb_t res_prime, ulong n, const arb_t x, slong prec)
67 {
68 arb_t xsub1, x2sub1;
69 double xx, xxsub1, cancellation_zero, cancellation_one;
70 double cost_zero, cost_one, cost_asymp;
71 double log2x, log2u, tolerance, asymp_error;
72 double yy, log2nsy, log2k, size;
73 slong wp;
74 slong d, k, K_zero, K_one, K_asymp;
75 int basecase_ok;
76
77 if (!arb_is_finite(x) || n > UWORD_MAX / 4)
78 {
79 if (res != NULL)
80 arb_indeterminate(res);
81 if (res_prime != NULL)
82 arb_indeterminate(res_prime);
83 return;
84 }
85
86 if (arf_sgn(arb_midref(x)) < 0)
87 {
88 arb_t t;
89 arb_init(t);
90 arb_neg(t, x);
91 arb_hypgeom_legendre_p_ui(res, res_prime, n, t, prec);
92 if (n % 2 == 1 && res != NULL)
93 arb_neg(res, res);
94 if (n % 2 == 0 && res_prime != NULL)
95 arb_neg(res_prime, res_prime);
96 arb_clear(t);
97 return;
98 }
99
100 if (arb_is_one(x) && n < UWORD_MAX)
101 {
102 if (res != NULL)
103 arb_set(res, x);
104 if (res_prime != NULL)
105 {
106 arb_set_ui(res_prime, n);
107 arb_mul_ui(res_prime, res_prime, n + 1, prec);
108 arb_mul_2exp_si(res_prime, res_prime, -1);
109 }
110 return;
111 }
112
113 if (n == 0)
114 {
115 if (res != NULL) arb_one(res);
116 if (res_prime != NULL) arb_zero(res_prime);
117 return;
118 }
119
120 if (n == 1)
121 {
122 if (res != NULL) arb_set_round(res, x, prec);
123 if (res_prime != NULL) arb_one(res_prime);
124 return;
125 }
126
127 xx = arf_get_d(arb_midref(x), ARF_RND_UP);
128
129 /* Use basecase recurrence? */
130 /* The following tests are not very elegant, and not completely accurate
131 either, but they are fast in the common case. */
132 if (res_prime != NULL)
133 {
134 basecase_ok = ((xx < 0.999999 && n < 10 && prec < 2000) ||
135 (xx < 0.999999 && n < 50 && prec < 1000) ||
136 (xx < 0.9999 && n < 100 && prec < 1000) ||
137 (xx < 0.999 && n < 350 && prec < 1000) ||
138 (xx < 0.9 && n < 400 && prec < 1000))
139 && ((xx > 0.00001 && n < 10 && prec < 2000) ||
140 (xx > 0.00001 && n < 60 && prec < 1000) ||
141 (xx > 0.01 && n < 200 && prec < 1000) ||
142 (xx > 0.1 && n < 400 && prec < 1000));
143
144 /* the recurrence also performs better when n ~= prec */
145 if (!basecase_ok)
146 basecase_ok = (xx > 0.1 && xx < 0.99 && n < 800
147 && prec > 0.4 * n && prec < 1.5 * n);
148 }
149 else if (prec < 500)
150 {
151 basecase_ok = ((xx < 0.999999 && n < 20) ||
152 (xx < 0.999 && n < 60) ||
153 (xx < 0.9 && n < 100))
154 && ((xx > 0.00001 && n < 20) ||
155 (xx > 0.01 && n < 60) ||
156 (xx > 0.1 && n < 100));
157
158 if (!basecase_ok)
159 basecase_ok = (xx > 0.1 && xx < 0.99 && n < 300
160 && prec > 0.4 * n && prec < 1.5 * n);
161 }
162 else
163 {
164 basecase_ok = 0;
165 }
166
167 if (basecase_ok)
168 {
169 mag_t t;
170 mag_init(t);
171 arb_get_mag(t, x);
172 if (mag_cmp_2exp_si(t, 0) >= 0)
173 basecase_ok = 0;
174 mag_clear(t);
175 }
176
177 if (basecase_ok)
178 {
179 arb_hypgeom_legendre_p_ui_rec(res, res_prime, n, x, prec);
180 return;
181 }
182
183 arb_init(xsub1);
184 arb_init(x2sub1);
185
186 arb_sub_ui(xsub1, x, 1, prec + 10);
187
188 arb_mul(x2sub1, x, x, 2 * prec);
189 arb_sub_ui(x2sub1, x2sub1, 1, prec + 10);
190 arb_neg(x2sub1, x2sub1);
191
192 /* use series at 1 unless |x| < 1-eps */
193 if (!arb_is_negative(xsub1) ||
194 arf_cmp_d(arb_midref(xsub1), ldexp(1.0, -2 * FLINT_BIT_COUNT(n))) > 0)
195 {
196 if (arf_cmp_d(arb_midref(xsub1), 2.0) >= 0)
197 {
198 if (n < 10000.0 * prec && n < UWORD_MAX / 4)
199 K_one = n + 1;
200 else
201 K_one = 1;
202 }
203 else /* check for early convergence */
204 {
205 xxsub1 = arf_get_d(arb_midref(xsub1), ARF_RND_UP);
206 log2u = log(fabs(xxsub1) * 0.5) * 1.44269504088896;
207 if (log2u < -30)
208 log2u = arf_abs_bound_lt_2exp_si(arb_midref(xsub1)) - 1.0;
209
210 K_one = n + 1;
211 K_one = FLINT_MIN(K_one, 100000.0 * prec);
212 K_one = FLINT_MIN(K_one, UWORD_MAX * 0.25);
213
214 size = 0.0;
215
216 if (n * (2.0 + log2u) < -prec)
217 {
218 for (k = 1; k < K_one; k = FLINT_MAX(k+1, k*1.05))
219 {
220 size = log2_bin_uiui_fast(n, k)
221 + log2_bin_uiui_fast(n + k, k) + k * log2u;
222
223 if (size < -prec)
224 {
225 K_one = k;
226 break;
227 }
228 }
229 }
230 }
231
232 arb_hypgeom_legendre_p_ui_one(res, res_prime, n, x, K_one, prec);
233 }
234 else /* guaranteed to have |x| < 1 */
235 {
236 cost_zero = 1e100;
237 cost_one = 1e100;
238 cost_asymp = 1e100;
239
240 xx = FLINT_MAX(xx, 1e-50);
241 xxsub1 = arf_get_d(arb_midref(xsub1), ARF_RND_UP);
242
243 /* Estimate cancellation for series expansion at 0. */
244 /* |P_n(xi)| ~= (x+sqrt(1+x^2))^n. */
245 cancellation_zero = n * log(xx + sqrt(1.0 + xx * xx)) * 1.44269504088896;
246 cancellation_zero = FLINT_MIN(cancellation_zero, 1.272 * n);
247 cancellation_zero = FLINT_MAX(cancellation_zero, 0.0);
248
249 /* Estimate cancellation for series expansion at 1. */
250 /* For x >= 1, P_n(x) ~= I_0(n*sqrt(2(x-1))) ~= exp(n*sqrt(2(x-1))) */
251 if (xxsub1 >= 0.0)
252 {
253 cancellation_one = 0.0;
254 }
255 else
256 {
257 cancellation_one = n * sqrt(2.0*fabs(xxsub1)) * 1.44269504088896;
258 cancellation_one = FLINT_MIN(cancellation_one, 2.0 * n);
259 cancellation_one = FLINT_MAX(cancellation_one, 0.0);
260 }
261
262 d = n / 2;
263 K_zero = d + 1;
264 K_one = n + 1;
265 K_asymp = 1;
266 asymp_error = 0.0;
267
268 wp = 1.01 * prec + FLINT_BIT_COUNT(n);
269 tolerance = -wp;
270
271 /* adjust for relative tolerance near 0 */
272 if (n % 2)
273 {
274 tolerance += arf_abs_bound_lt_2exp_si(arb_midref(x));
275 }
276
277 if (n > 10)
278 {
279 /* look for early truncation of series at 1 */
280 log2u = log(fabs(xxsub1) * 0.5) * 1.44269504088896;
281 if (log2u < -30)
282 log2u = arf_abs_bound_lt_2exp_si(arb_midref(xsub1)) - 1.0;
283
284 log2x = log(fabs(xx)) * 1.44269504088896;
285 if (log2x < -30)
286 log2x = arf_abs_bound_lt_2exp_si(arb_midref(x));
287
288 if (n * (2.0 + log2u) < tolerance)
289 {
290 for (k = 1; k < K_one; k = FLINT_MAX(k+1, k*1.05))
291 {
292 size = log2_bin_uiui_fast(n, k)
293 + log2_bin_uiui_fast(n + k, k) + k * log2u;
294
295 if (size < tolerance)
296 {
297 K_one = k;
298 break;
299 }
300 }
301 }
302
303 /* look for early truncation of series at 0 */
304 if (n * (1.0 + log2x) < tolerance)
305 {
306 for (k = 1; k < K_zero; k = FLINT_MAX(k+1, k*1.05))
307 {
308 size = log2_bin_uiui_fast(n, d - k)
309 + log2_bin_uiui_fast(n+1+2*k, n) - n + 2.0 * k * log2x;
310
311 if (size < tolerance)
312 {
313 K_zero = k;
314 break;
315 }
316 }
317 }
318
319 /* look for possible convergence of asymptotic series */
320 /* requires information about y = sqrt(1-x^2) */
321 yy = arf_get_d(arb_midref(x2sub1), ARF_RND_DOWN);
322 yy = sqrt(FLINT_MAX(0.0, yy));
323 log2nsy = log(2.0 * n * yy) * 1.44269504088896;
324
325 cost_zero = (prec + cancellation_zero) * K_zero;
326 cost_one = (prec + cancellation_one) * K_one;
327
328 for (k = 1; k < n &&
329 prec * k < FLINT_MIN(cost_zero, cost_one);
330 k = FLINT_MAX(k + 1, k * 1.05))
331 {
332 /* todo: better account for prefactor in the asymptotic series? */
333 log2k = log(k) * 1.44269504088896;
334 size = 3.0 + k * (log2k - 1.43); /* estimate K! */
335 size -= k * log2nsy; /* 1/(2n sin(theta))^K */
336
337 if (size < asymp_error)
338 {
339 asymp_error = size;
340 K_asymp = k;
341 }
342
343 if (size < tolerance)
344 {
345 break;
346 }
347 }
348 }
349
350 cost_zero = (prec + cancellation_zero) * K_zero;
351 cost_one = (prec + cancellation_one) * K_one;
352 cost_asymp = (prec + 0.0) * K_asymp * 2.0;
353
354 #if 0
355 printf("zero: K = %ld, cost = %g, cancel %g\n", K_zero, cost_zero, cancellation_zero);
356 printf("one: K = %ld, cost = %g, cancel %g\n", K_one, cost_one, cancellation_one);
357 printf("asymp: K = %ld, cost = %g, error = %f (tol = %f)\n", K_asymp, cost_asymp, asymp_error, tolerance);
358 #endif
359
360 if (asymp_error < tolerance && cost_asymp < FLINT_MIN(cost_zero, cost_one))
361 {
362 arb_hypgeom_legendre_p_ui_asymp(res, res_prime, n, x, K_asymp, wp);
363 }
364 else if (FLINT_MIN(cost_zero, cost_one) < (1e6 * prec) * prec && n < UWORD_MAX / 4)
365 {
366 mag_t err1, err2, xrad;
367 arb_t xmid;
368
369 mag_init(err1);
370 mag_init(err2);
371 mag_init(xrad);
372 arb_init(xmid);
373
374 arf_set(arb_midref(xmid), arb_midref(x));
375 mag_zero(arb_radref(xmid));
376 mag_set(xrad, arb_radref(x));
377
378 arb_hypgeom_legendre_p_ui_deriv_bound(err1, err2, n, x, x2sub1);
379
380 if (cost_zero < cost_one)
381 arb_hypgeom_legendre_p_ui_zero(res, res_prime, n, xmid, K_zero, wp + cancellation_zero);
382 else
383 arb_hypgeom_legendre_p_ui_one(res, res_prime, n, xmid, K_one, wp + cancellation_one);
384
385 if (res != NULL)
386 {
387 mag_mul(err1, err1, xrad);
388 arb_add_error_mag(res, err1);
389 arb_set_round(res, res, prec);
390 }
391
392 if (res_prime != NULL)
393 {
394 mag_mul(err2, err2, xrad);
395 arb_add_error_mag(res_prime, err2);
396 arb_set_round(res_prime, res_prime, prec);
397 }
398
399 mag_clear(err1);
400 mag_clear(err2);
401 mag_clear(xrad);
402 arb_clear(xmid);
403 }
404 else if (asymp_error < -2.0)
405 {
406 /* todo -- clamp to [-1,1]? */
407 arb_hypgeom_legendre_p_ui_asymp(res, res_prime, n, x, K_asymp, wp);
408 }
409 else
410 {
411 if (res != NULL)
412 {
413 arf_zero(arb_midref(res));
414 mag_one(arb_radref(res));
415 }
416
417 if (res_prime != NULL)
418 {
419 arf_zero(arb_midref(res_prime));
420 mag_set_ui(arb_radref(res_prime), n);
421 mag_add_ui(arb_radref(res_prime), arb_radref(res_prime), 1);
422 mag_mul_ui(arb_radref(res_prime), arb_radref(res_prime), n);
423 mag_mul_2exp_si(arb_radref(res_prime), arb_radref(res_prime), -1);
424 }
425 }
426 }
427
428 arb_clear(xsub1);
429 arb_clear(x2sub1);
430 }
431
432