1 /*
2 Copyright (C) 2019 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 /* need special case for integer l and z = 0 since the recurrence relations break down */
15 static void
_acb_hypgeom_coulomb_f_int_jet(acb_ptr F,const acb_t l,const acb_t eta,const acb_t z,slong len,slong prec)16 _acb_hypgeom_coulomb_f_int_jet(acb_ptr F, const acb_t l, const acb_t eta, const acb_t z, slong len, slong prec)
17 {
18 acb_poly_struct a[1];
19 acb_poly_struct b[2];
20 acb_poly_t zx, M, zxi;
21 acb_t t;
22 slong k;
23 int real;
24
25 if (arf_cmp_si(arb_midref(acb_realref(l)), -1) < 0)
26 {
27 _acb_vec_indeterminate(F, len);
28 return;
29 }
30
31 /* http://fungrim.org/entry/2a2f18/ */
32 /* F = C * (z+x)^(l+1) e^(-+ i (z+x)) M(l + 1 -+ i eta, 2l+2, +- 2 i (z+x)) */
33
34 acb_poly_init(a);
35 acb_poly_init(b);
36 acb_poly_init(b + 1);
37 acb_poly_init(zx);
38 acb_poly_init(M);
39 acb_poly_init(zxi);
40 acb_init(t);
41
42 acb_poly_set_coeff_acb(zx, 0, z);
43 acb_poly_set_coeff_si(zx, 1, 1);
44
45 acb_div_onei(t, eta);
46 acb_add(t, t, l, prec);
47 acb_add_ui(t, t, 1, prec);
48 acb_poly_set_acb(a, t);
49
50 acb_add_ui(t, l, 1, prec);
51 acb_mul_2exp_si(t, t, 1);
52 acb_poly_set_acb(b, t);
53
54 acb_poly_one(b + 1);
55
56 acb_onei(t);
57 acb_mul_2exp_si(t, t, 1);
58 acb_poly_scalar_mul(zxi, zx, t, prec);
59
60 acb_hypgeom_pfq_series_direct(M, a, 1, b, 2, zxi, 1, -1, len, prec);
61
62 acb_poly_scalar_mul_2exp_si(zxi, zxi, -1);
63 acb_poly_neg(zxi, zxi);
64 acb_poly_exp_series(zxi, zxi, len, prec);
65 acb_poly_mullow(M, M, zxi, len, prec);
66
67 acb_add_ui(t, l, 1, prec);
68 acb_poly_pow_acb_series(zxi, zx, t, len, prec);
69
70 acb_poly_mullow(M, M, zxi, len, prec);
71
72 {
73 /* C = 2^l exp((-pi eta + lu + lv)/2) */
74 acb_t C, lu, lv;
75
76 acb_init(C);
77 acb_init(lu);
78 acb_init(lv);
79
80 acb_add_ui(lu, l, 1, prec);
81 acb_mul_onei(t, eta);
82 acb_add(lu, lu, t, prec);
83
84 acb_add_ui(lv, l, 1, prec);
85 acb_div_onei(t, eta);
86 acb_add(lv, lv, t, prec);
87
88 acb_lgamma(lu, lu, prec);
89 acb_lgamma(lv, lv, prec);
90
91 acb_const_pi(t, prec);
92 acb_add(C, lu, lv, prec);
93 acb_submul(C, t, eta, prec);
94 acb_mul_2exp_si(C, C, -1);
95 acb_exp(C, C, prec);
96
97 acb_set_ui(t, 2);
98 acb_pow(t, t, l, prec);
99 acb_mul(C, C, t, prec);
100
101 acb_poly_scalar_mul(M, M, C, prec);
102
103 acb_clear(C);
104 acb_clear(lu);
105 acb_clear(lv);
106 }
107
108 real = acb_is_real(z) && acb_is_real(eta);
109
110 for (k = 0; k < len; k++)
111 {
112 acb_poly_get_coeff_acb(F + k, M, k);
113 if (real)
114 arb_zero(acb_imagref(F + k));
115 }
116
117 acb_poly_clear(a);
118 acb_poly_clear(b);
119 acb_poly_clear(b + 1);
120 acb_poly_clear(zx);
121 acb_poly_clear(M);
122 acb_poly_clear(zxi);
123 acb_clear(t);
124 }
125
126 static void
_acb_hypgeom_coulomb_jet(acb_ptr F,acb_ptr G,acb_ptr Hpos,acb_ptr Hneg,const acb_t l,const acb_t eta,const acb_t z,slong len,slong prec)127 _acb_hypgeom_coulomb_jet(acb_ptr F, acb_ptr G, acb_ptr Hpos, acb_ptr Hneg, const acb_t l, const acb_t eta, const acb_t z, slong len, slong prec)
128 {
129 acb_t l1, t, R, S;
130
131 if (len <= 0)
132 return;
133
134 if (len == 1)
135 {
136 acb_hypgeom_coulomb(F, G, Hpos, Hneg, l, eta, z, prec);
137 return;
138 }
139
140 if (acb_contains_zero(z))
141 {
142 if (F != NULL)
143 {
144 if (acb_is_int(l))
145 _acb_hypgeom_coulomb_f_int_jet(F, l, eta, z, len, prec);
146 else
147 _acb_vec_indeterminate(F, len);
148 }
149
150 if (G != NULL) _acb_vec_indeterminate(G, len);
151 if (Hpos != NULL) _acb_vec_indeterminate(Hpos, len);
152 if (Hneg != NULL) _acb_vec_indeterminate(Hneg, len);
153 return;
154 }
155
156
157 acb_init(l1);
158 acb_init(t);
159 acb_init(R);
160 acb_init(S);
161
162 acb_add_ui(l1, l, 1, prec);
163
164 acb_hypgeom_coulomb(F, G, Hpos, Hneg, l, eta, z, prec);
165
166 /* todo: somehow recycle the gamma function values for the two evaluations? */
167 acb_hypgeom_coulomb((F == NULL) ? NULL : (F + 1),
168 (G == NULL) ? NULL : (G + 1),
169 (Hpos == NULL) ? NULL : (Hpos + 1),
170 (Hneg == NULL) ? NULL : (Hneg + 1),
171 l1, eta, z, prec);
172
173 /* First derivatives:
174 http://fungrim.org/entry/a51a4b/, http://fungrim.org/entry/2fec14/ */
175
176 /* R_l = (2l+1) C_l / C_{l+1} = sqrt(l+i eta) sqrt(l-i eta) / l */
177 if (acb_is_real(l) && acb_is_real(eta) && arb_is_nonzero(acb_realref(eta)))
178 {
179 acb_mul(R, l1, l1, prec);
180 acb_addmul(R, eta, eta, prec);
181 acb_sqrt(R, R, prec);
182 }
183 else
184 {
185 acb_mul_onei(t, eta);
186 acb_add(t, t, l1, prec);
187 acb_sqrt(t, t, prec);
188 acb_div_onei(R, eta);
189 acb_add(R, R, l1, prec);
190 acb_sqrt(R, R, prec);
191 acb_mul(R, R, t, prec);
192 }
193 acb_div(R, R, l1, prec);
194
195 acb_div(S, l1, z, prec);
196 acb_div(t, eta, l1, prec);
197 acb_add(S, S, t, prec);
198
199 /* todo: fix regular F at origin */
200 if (F != NULL)
201 {
202 acb_mul(F + 1, F + 1, R, prec);
203 acb_neg(F + 1, F + 1);
204 acb_addmul(F + 1, F, S, prec);
205 }
206
207 if (G != NULL)
208 {
209 acb_mul(G + 1, G + 1, R, prec);
210 acb_neg(G + 1, G + 1);
211 acb_addmul(G + 1, G, S, prec);
212 }
213
214 if (Hpos != NULL)
215 {
216 acb_mul(Hpos + 1, Hpos + 1, R, prec);
217 acb_neg(Hpos + 1, Hpos + 1);
218 acb_addmul(Hpos + 1, Hpos, S, prec);
219 }
220
221 if (Hneg != NULL)
222 {
223 acb_mul(Hneg + 1, Hneg + 1, R, prec);
224 acb_neg(Hneg + 1, Hneg + 1);
225 acb_addmul(Hneg + 1, Hneg, S, prec);
226 }
227
228 if (len >= 3)
229 {
230 acb_t q, q2, w, w2;
231
232 acb_init(q);
233 acb_init(q2);
234 acb_init(w);
235 acb_init(w2);
236
237 acb_inv(w, z, prec);
238 acb_mul(w2, w, w, prec);
239
240 /* http://fungrim.org/entry/07a654/ */
241 /* F''/2 = q F, q = (2eta/z + l(l+1)/z^2 - 1)/2 */
242
243 acb_mul(q, l, l1, prec);
244 acb_mul(q, q, w2, prec);
245 acb_mul_2exp_si(q2, eta, 1);
246 acb_addmul(q, q2, w, prec);
247 acb_sub_ui(q, q, 1, prec);
248 acb_mul_2exp_si(q, q, -1);
249
250 if (F != NULL) acb_mul(F + 2, F, q, prec);
251 if (G != NULL) acb_mul(G + 2, G, q, prec);
252 if (Hpos != NULL) acb_mul(Hpos + 2, Hpos, q, prec);
253 if (Hneg != NULL) acb_mul(Hneg + 2, Hneg, q, prec);
254
255 /* http://fungrim.org/entry/faa118/ */
256 /* F'''/6 = (2qF' - q2 F)/6, q2 = 2(eta + l(l+1)/z)/z^2 */
257 if (len >= 4)
258 {
259 acb_mul_2exp_si(q, q, 1);
260 acb_mul(q2, l, l1, prec);
261 acb_mul(q2, q2, w, prec);
262 acb_add(q2, q2, eta, prec);
263 acb_mul_2exp_si(q2, q2, 1);
264 acb_mul(q2, q2, w2, prec);
265
266 if (F != NULL)
267 {
268 acb_mul(F + 3, F + 1, q, prec);
269 acb_submul(F + 3, F + 0, q2, prec);
270 acb_div_ui(F + 3, F + 3, 6, prec);
271 }
272
273 if (G != NULL)
274 {
275 acb_mul(G + 3, G + 1, q, prec);
276 acb_submul(G + 3, G, q2, prec);
277 acb_div_ui(G + 3, G + 3, 6, prec);
278 }
279
280 if (Hpos != NULL)
281 {
282 acb_mul(Hpos + 3, Hpos + 1, q, prec);
283 acb_submul(Hpos + 3, Hpos, q2, prec);
284 acb_div_ui(Hpos + 3, Hpos + 3, 6, prec);
285 }
286
287 if (Hneg != NULL)
288 {
289 acb_mul(Hneg + 3, Hneg + 1, q, prec);
290 acb_submul(Hneg + 3, Hneg, q2, prec);
291 acb_div_ui(Hneg + 3, Hneg + 3, 6, prec);
292 }
293 }
294
295 /* http://fungrim.org/entry/eca10b/ */
296 if (len >= 5)
297 {
298 slong k;
299 acb_ptr s;
300 s = _acb_vec_init(4);
301
302 /*
303 s4 = -(k(k+7) + 12)
304 s3 = 2 (k(k+5) + 6) z
305 s2 = (k(k+3) + z^2 - 2z*eta - l(l+1) + 2)
306 s1 = 2(z-eta)
307 s0 = 1
308 F[k+4] = (s0*F[k] + s1*F[k+1] + s2*F[k+2] + s3*F[k+3]) / s4 / z^2
309 */
310 acb_sub(s + 1, z, eta, prec);
311 acb_mul_2exp_si(s + 1, s + 1, 1);
312
313 acb_mul(q, z, z, prec);
314 acb_mul(q2, eta, z, prec);
315 acb_mul_2exp_si(q2, q2, 1);
316 acb_sub(q, q, q2, prec);
317 acb_submul(q, l, l1, prec);
318
319 for (k = 0; k + 4 < len; k++)
320 {
321 acb_mul_si(s + 3, z, 2 * (k * (k + 5) + 6), prec);
322 acb_add_si(s + 2, q, k * (k + 3) + 2, prec);
323
324 if (F != NULL)
325 {
326 acb_dot(F + k + 4, F + k, 0, F + k + 1, 1, s + 1, 1, 3, prec);
327 acb_div_si(F + k + 4, F + k + 4, -(k * (k + 7) + 12), prec);
328 acb_mul(F + k + 4, F + k + 4, w2, prec);
329 }
330
331 if (G != NULL)
332 {
333 acb_dot(G + k + 4, G + k, 0, G + k + 1, 1, s + 1, 1, 3, prec);
334 acb_div_si(G + k + 4, G + k + 4, -(k * (k + 7) + 12), prec);
335 acb_mul(G + k + 4, G + k + 4, w2, prec);
336 }
337
338 if (Hpos != NULL)
339 {
340 acb_dot(Hpos + k + 4, Hpos + k, 0, Hpos + k + 1, 1, s + 1, 1, 3, prec);
341 acb_div_si(Hpos + k + 4, Hpos + k + 4, -(k * (k + 7) + 12), prec);
342 acb_mul(Hpos + k + 4, Hpos + k + 4, w2, prec);
343 }
344
345 if (Hneg != NULL)
346 {
347 acb_dot(Hneg + k + 4, Hneg + k, 0, Hneg + k + 1, 1, s + 1, 1, 3, prec);
348 acb_div_si(Hneg + k + 4, Hneg + k + 4, -(k * (k + 7) + 12), prec);
349 acb_mul(Hneg + k + 4, Hneg + k + 4, w2, prec);
350 }
351 }
352
353 _acb_vec_clear(s, 4);
354 }
355
356 acb_clear(q);
357 acb_clear(q2);
358 acb_clear(w);
359 acb_clear(w2);
360 }
361
362 acb_clear(l1);
363 acb_clear(t);
364 acb_clear(R);
365 acb_clear(S);
366 }
367
368 void
acb_hypgeom_coulomb_jet(acb_ptr F,acb_ptr G,acb_ptr Hpos,acb_ptr Hneg,const acb_t l,const acb_t eta,const acb_t z,slong len,slong prec)369 acb_hypgeom_coulomb_jet(acb_ptr F, acb_ptr G, acb_ptr Hpos, acb_ptr Hneg, const acb_t l, const acb_t eta, const acb_t z, slong len, slong prec)
370 {
371 acb_t t; /* to allow aliasing with z */
372 acb_init(t);
373 acb_set(t, z);
374 _acb_hypgeom_coulomb_jet(F, G, Hpos, Hneg, l, eta, t, len, prec);
375 acb_clear(t);
376 }
377
378