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 void
acb_hypgeom_bessel_k_asymp(acb_t res,const acb_t nu,const acb_t z,int scaled,slong prec)15 acb_hypgeom_bessel_k_asymp(acb_t res, const acb_t nu, const acb_t z, int scaled, slong prec)
16 {
17 acb_t t, a, b, w;
18
19 acb_init(t);
20 acb_init(a);
21 acb_init(b);
22 acb_init(w);
23
24 acb_one(a);
25 acb_mul_2exp_si(a, a, -1);
26 acb_add(a, a, nu, prec);
27
28 acb_mul_2exp_si(b, nu, 1);
29 acb_add_ui(b, b, 1, prec);
30
31 acb_mul_2exp_si(w, z, 1);
32
33 acb_hypgeom_u_asymp(t, a, b, w, -1, prec);
34
35 if (!scaled)
36 {
37 acb_neg(a, z);
38 acb_exp(a, a, prec);
39 acb_mul(t, t, a, prec);
40 }
41
42 acb_mul_2exp_si(w, z, 1);
43 acb_rsqrt(w, w, prec);
44 acb_mul(res, t, w, prec);
45
46 arb_const_sqrt_pi(acb_realref(w), prec);
47 acb_mul_arb(res, res, acb_realref(w), prec);
48
49 acb_clear(t);
50 acb_clear(a);
51 acb_clear(b);
52 acb_clear(w);
53 }
54
55 void
acb_hypgeom_bessel_k_0f1_series(acb_poly_t res,const acb_poly_t nu,const acb_poly_t z,int scaled,slong len,slong prec)56 acb_hypgeom_bessel_k_0f1_series(acb_poly_t res,
57 const acb_poly_t nu, const acb_poly_t z,
58 int scaled, slong len, slong prec)
59 {
60 acb_poly_t s, u, A, B;
61 acb_poly_struct b[2];
62 arb_t c;
63 slong wlen;
64 int singular;
65
66 acb_poly_init(s);
67 acb_poly_init(u);
68 acb_poly_init(A);
69 acb_poly_init(B);
70 acb_poly_init(b + 0);
71 acb_poly_init(b + 1);
72 arb_init(c);
73
74 singular = (nu->length == 0) || acb_is_int(nu->coeffs);
75 wlen = len + (singular != 0);
76
77 /* A = (z/2)^nu, B = 1/A */
78 acb_poly_scalar_mul_2exp_si(A, z, -1);
79 acb_poly_pow_series(A, A, nu, wlen, prec);
80 acb_poly_inv_series(B, A, wlen, prec);
81
82 acb_poly_mullow(u, z, z, wlen, prec);
83 acb_poly_scalar_mul_2exp_si(u, u, -2);
84
85 acb_poly_one(b + 1);
86 acb_poly_add_si(b + 0, nu, 1, prec);
87 acb_hypgeom_pfq_series_direct(s, NULL, 0, b, 2, u, 1, -1, wlen, prec);
88 acb_poly_mullow(A, A, s, wlen, prec);
89
90 acb_poly_add_si(b + 0, nu, -1, prec);
91 acb_poly_neg(b + 0, b + 0);
92 acb_hypgeom_pfq_series_direct(s, NULL, 0, b, 2, u, 1, -1, wlen, prec);
93 acb_poly_mullow(B, B, s, wlen, prec);
94
95 acb_poly_sub(A, B, A, prec);
96 acb_poly_scalar_mul_2exp_si(A, A, -1);
97
98 /* multiply by pi csc(pi nu) */
99 acb_poly_sin_pi_series(B, nu, wlen, prec);
100
101 if (singular)
102 {
103 acb_poly_shift_right(A, A, 1);
104 acb_poly_shift_right(B, B, 1);
105 }
106
107 if (scaled)
108 {
109 acb_poly_exp_series(u, z, len, prec);
110 acb_poly_mullow(A, A, u, len, prec);
111 }
112
113 acb_poly_div_series(res, A, B, len, prec);
114
115 arb_const_pi(c, prec);
116 _acb_vec_scalar_mul_arb(res->coeffs, res->coeffs, res->length, c, prec);
117
118 acb_poly_clear(s);
119 acb_poly_clear(u);
120 acb_poly_clear(A);
121 acb_poly_clear(B);
122 acb_poly_clear(b + 0);
123 acb_poly_clear(b + 1);
124 arb_clear(c);
125 }
126
127 void
acb_hypgeom_bessel_k_0f1(acb_t res,const acb_t nu,const acb_t z,int scaled,slong prec)128 acb_hypgeom_bessel_k_0f1(acb_t res, const acb_t nu, const acb_t z, int scaled, slong prec)
129 {
130 if (acb_is_int(nu))
131 {
132 acb_poly_t nux, zx, rx;
133
134 acb_poly_init(nux);
135 acb_poly_init(zx);
136 acb_poly_init(rx);
137
138 acb_poly_set_coeff_acb(nux, 0, nu);
139 acb_poly_set_coeff_si(nux, 1, 1);
140 acb_poly_set_acb(zx, z);
141
142 acb_hypgeom_bessel_k_0f1_series(rx, nux, zx, scaled, 1, prec);
143
144 acb_poly_get_coeff_acb(res, rx, 0);
145
146 acb_poly_clear(nux);
147 acb_poly_clear(zx);
148 acb_poly_clear(rx);
149 }
150 else
151 {
152 acb_t t, u, v, w;
153 acb_struct b[2];
154
155 acb_init(t);
156 acb_init(u);
157 acb_init(v);
158 acb_init(w);
159 acb_init(b + 0);
160 acb_init(b + 1);
161
162 /* u = 0F1(1+nu), v = 0F1(1-nu) */
163 acb_mul(t, z, z, prec);
164 acb_mul_2exp_si(t, t, -2);
165 acb_add_ui(b, nu, 1, prec);
166 acb_one(b + 1);
167 acb_hypgeom_pfq_direct(u, NULL, 0, b, 2, t, -1, prec);
168 acb_sub_ui(b, nu, 1, prec);
169 acb_neg(b, b);
170 acb_hypgeom_pfq_direct(v, NULL, 0, b, 2, t, -1, prec);
171
172 /* v = v * gamma(nu) / (z/2)^nu */
173 acb_mul_2exp_si(t, z, -1);
174 acb_pow(t, t, nu, prec);
175 acb_gamma(w, nu, prec);
176 acb_mul(v, v, w, prec);
177 acb_div(v, v, t, prec);
178
179 /* u = u * t * pi / (gamma(nu) * nu * sin(pi nu)) */
180 acb_mul(u, u, t, prec);
181 acb_const_pi(t, prec);
182 acb_mul(u, u, t, prec);
183 acb_sin_pi(t, nu, prec);
184 acb_mul(t, t, w, prec);
185 acb_mul(t, t, nu, prec);
186 acb_div(u, u, t, prec);
187
188 acb_sub(v, v, u, prec);
189 acb_mul_2exp_si(v, v, -1);
190
191 if (scaled)
192 {
193 acb_exp(t, z, prec);
194 acb_mul(v, v, t, prec);
195 }
196
197 acb_set(res, v);
198
199 acb_clear(t);
200 acb_clear(u);
201 acb_clear(v);
202 acb_clear(w);
203 acb_clear(b + 0);
204 acb_clear(b + 1);
205 }
206 }
207
208 void
acb_hypgeom_bessel_k(acb_t res,const acb_t nu,const acb_t z,slong prec)209 acb_hypgeom_bessel_k(acb_t res, const acb_t nu, const acb_t z, slong prec)
210 {
211 mag_t zmag;
212
213 mag_init(zmag);
214 acb_get_mag(zmag, z);
215
216 if (mag_cmp_2exp_si(zmag, 4) < 0 ||
217 (mag_cmp_2exp_si(zmag, 64) < 0 && 2 * mag_get_d(zmag) < prec))
218 acb_hypgeom_bessel_k_0f1(res, nu, z, 0, prec);
219 else
220 acb_hypgeom_bessel_k_asymp(res, nu, z, 0, prec);
221
222 mag_clear(zmag);
223 }
224
225 void
acb_hypgeom_bessel_k_scaled(acb_t res,const acb_t nu,const acb_t z,slong prec)226 acb_hypgeom_bessel_k_scaled(acb_t res, const acb_t nu, const acb_t z, slong prec)
227 {
228 mag_t zmag;
229
230 mag_init(zmag);
231 acb_get_mag(zmag, z);
232
233 if (mag_cmp_2exp_si(zmag, 4) < 0 ||
234 (mag_cmp_2exp_si(zmag, 64) < 0 && 2 * mag_get_d(zmag) < prec))
235 acb_hypgeom_bessel_k_0f1(res, nu, z, 1, prec);
236 else
237 acb_hypgeom_bessel_k_asymp(res, nu, z, 1, prec);
238
239 mag_clear(zmag);
240 }
241
242