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 #define UNROLL 4
15
16 static void
sum_rs_inner(arb_t s,arb_srcptr xpow,slong m,ulong n,slong K,slong prec)17 sum_rs_inner(arb_t s, arb_srcptr xpow, slong m, ulong n, slong K, slong prec)
18 {
19 slong j, k, khi, klo, u, r;
20 ulong d;
21 slong sigma;
22 fmpz * c;
23
24 sigma = n % 2 ? 1 : -1;
25 d = n / 2;
26
27 arb_zero(s);
28 c = _fmpz_vec_init(UNROLL + 1);
29
30 k = K - 1;
31 while (k >= 1)
32 {
33 u = FLINT_MIN(UNROLL, k);
34
35 khi = k;
36 klo = k - u + 1;
37
38 for (j = klo; j <= khi; j++)
39 {
40 ulong aa = (d - j + 1);
41 ulong bb = (2 * d + 2 * j + sigma);
42
43 if (j == klo)
44 fmpz_ui_mul_ui(c + khi - j, aa, bb);
45 else
46 fmpz_mul2_uiui(c + khi - j, c + khi - j + 1, aa, bb);
47 }
48
49 for (j = khi; j >= klo; j--)
50 {
51 ulong aa = (j);
52 ulong bb = (2 * j + sigma);
53
54 if (j == khi)
55 {
56 fmpz_ui_mul_ui(c + u, aa, bb);
57 }
58 else
59 {
60 fmpz_mul(c + khi - j, c + khi - j, c + u);
61 fmpz_mul2_uiui(c + u, c + u, aa, bb);
62 }
63 }
64
65 while (k >= klo)
66 {
67 r = k % m;
68 if (k == khi)
69 {
70 arb_add(s, s, xpow + r, prec);
71 arb_mul_fmpz(s, s, c + khi - k, prec);
72 }
73 else if (r == 0)
74 arb_add_fmpz(s, s, c + khi - k, prec);
75 else
76 arb_addmul_fmpz(s, xpow + r, c + khi - k, prec);
77
78 if (r == 0 && k != 0)
79 arb_mul(s, s, xpow + m, prec);
80
81 k--;
82 }
83
84 arb_div_fmpz(s, s, c + u, prec);
85 }
86
87 _fmpz_vec_clear(c, UNROLL + 1);
88 }
89
90 static void
_arb_hypgeom_legendre_p_ui_zero(arb_t res,ulong n,const arb_t x,arb_srcptr xpow,slong m,slong K,slong prec)91 _arb_hypgeom_legendre_p_ui_zero(arb_t res, ulong n,
92 const arb_t x, arb_srcptr xpow, slong m, slong K, slong prec)
93 {
94 arb_t s;
95 slong d, prec2;
96 mag_t u, a, t, err;
97
98 d = n / 2;
99
100 arb_init(s);
101 mag_init(u);
102 mag_init(a);
103 mag_init(t);
104 mag_init(err);
105
106 K = FLINT_MIN(K, d + 1);
107 sum_rs_inner(s, xpow, m, n, K, prec);
108
109 prec2 = arb_rel_accuracy_bits(s);
110 if (prec2 > prec)
111 prec2 = prec;
112 else
113 prec2 = FLINT_MAX(0, prec2) + 20;
114
115 arb_add_ui(s, s, 1, prec2);
116 if (n % 2 == 1)
117 arb_mul(s, s, x, prec2);
118 arb_swap(res, s);
119
120 if (d % 2 == 1)
121 arb_neg(res, res);
122
123 if (n % 2 == 0)
124 {
125 arb_hypgeom_central_bin_ui(s, d, prec2);
126 arb_mul(res, res, s, prec2);
127 arb_mul_2exp_si(res, res, -n);
128 }
129 else
130 {
131 arb_hypgeom_central_bin_ui(s, d + 1, prec2);
132 arb_mul(res, res, s, prec2);
133 arb_mul_ui(res, res, d + 1, prec2);
134 arb_mul_2exp_si(res, res, -n);
135 }
136
137 if (K < d + 1)
138 {
139 mag_bin_uiui(err, n, d - K);
140 mag_bin_uiui(t, n + 2 * K + (n % 2), n);
141 mag_mul(err, err, t);
142 arb_get_mag(t, x);
143 mag_pow_ui(t, t, 2 * K + (n % 2));
144 mag_mul(err, err, t);
145 mag_mul_2exp_si(err, err, -n);
146
147 arb_get_mag(t, x);
148 mag_mul(t, t, t);
149 mag_mul_ui(t, t, d - K + 1);
150 mag_mul_ui(t, t, 2 * d + 2 * K + (n % 2 ? 1 : -1));
151 mag_div_ui(t, t, K);
152 mag_div_ui(t, t, 2 * K + (n % 2 ? 1 : -1));
153 mag_geom_series(t, t, 0);
154 mag_mul(err, err, t);
155
156 arb_add_error_mag(res, err);
157 }
158
159 arb_clear(s);
160 mag_clear(u);
161 mag_clear(a);
162 mag_clear(t);
163 mag_clear(err);
164 }
165
166
167 void
arb_hypgeom_legendre_p_ui_zero(arb_t res,arb_t res2,ulong n,const arb_t x,slong K,slong prec)168 arb_hypgeom_legendre_p_ui_zero(arb_t res, arb_t res2, ulong n,
169 const arb_t x, slong K, slong prec)
170 {
171 arb_ptr xpow;
172 arb_t t, u, v;
173 slong m;
174
175 if (n == 0)
176 {
177 if (res != NULL) arb_one(res);
178 if (res2 != NULL) arb_zero(res2);
179 return;
180 }
181
182 /* overflow protection */
183 if (n > UWORD_MAX / 4)
184 {
185 if (res != NULL) arb_indeterminate(res);
186 if (res2 != NULL) arb_indeterminate(res2);
187 }
188
189 if (res == NULL)
190 {
191 arb_init(v);
192 arb_hypgeom_legendre_p_ui_zero(v, res2, n, x, K, prec);
193 arb_clear(v);
194 return;
195 }
196
197 arb_init(t);
198 arb_init(u);
199 arb_init(v);
200
201 K = FLINT_MIN(K, n / 2 + 1);
202
203 if (res2 != NULL)
204 m = n_sqrt(2 * K);
205 else
206 m = n_sqrt(K);
207
208 xpow = _arb_vec_init(m + 1);
209
210 arb_mul(v, x, x, prec);
211 arb_neg(v, v);
212 _arb_vec_set_powers(xpow, v, m + 1, prec);
213
214 /* todo: recycle prefactor */
215
216 if (res2 == NULL)
217 {
218 _arb_hypgeom_legendre_p_ui_zero(t, n, x, xpow, m, K, prec);
219 arb_set(res, t);
220 }
221 else
222 {
223 _arb_hypgeom_legendre_p_ui_zero(t, n, x, xpow, m, K, prec);
224 _arb_hypgeom_legendre_p_ui_zero(u, n - 1, x, xpow, m, K, prec);
225
226 /* P' = n (P[n-1] - x P) / (1 - x^2) */
227 arb_submul(u, t, x, prec);
228 arb_add_ui(v, v, 1, prec);
229 arb_div(u, u, v, prec);
230 arb_mul_ui(res2, u, n, prec);
231 arb_set(res, t);
232 }
233
234 _arb_vec_clear(xpow, m + 1);
235 arb_clear(t);
236 arb_clear(u);
237 arb_clear(v);
238 }
239
240