1 /*
2 Copyright (C) 2018 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 /*
15 https://dlmf.nist.gov/9.9
16
17 a_k ~ -T(3/8 pi (4k-1))
18 a'_k ~ -U(3/8 pi (4k-3))
19 b_k ~ -T(3/8 pi (4k-3))
20 b'_k ~ -U(3/8 pi (4k-1))
21
22 For a_k and b_k, the u^8 and u^10 truncations are known to give lower
23 bounds. [G. Pittaluga and L. Sacripante (1991) Inequalities for the
24 zeros of the Airy functions. SIAM J. Math. Anal. 22 (1), pp. 260–267.]
25
26 We don't have proofs for a'_k and b'_k. However, in that case, we can just
27 do a single interval Newton step to verify that we have isolated a
28 zero (the enclosure must be for the correct zero due to sandwiching).
29 */
30
31 #define AI 0
32 #define BI 1
33 #define AI_PRIME 2
34 #define BI_PRIME 3
35
36 static const double initial[4][10] = {{
37 -658118728906175.0,-1150655474581104.0,-1553899449042978.0,-1910288501594969.0,
38 -2236074816421182.0,-2539650438812533.0,-2826057838960988.0,
39 -3098624122012011.0,-3359689702679955.0,-3610979637739094.0},{
40 -330370902027041.0,-920730911234245.0,-1359731821477101.0,-1736658984124319.0,
41 -2076373934490092.0,-2390271103799312.0,-2684763040788193.0,
42 -2963907159065113.0,-3230475233555475.0,-3486466475611047.0},{
43 -286764727967452.0,-914286338795679.0,-1356737313209586.0,-1734816794389239.0,
44 -2075083421171399.0,-2389296605766914.0,-2683990299959380.0,
45 -2963272965051282.0,-3229941298662311.0,-3486008018531685.0},{
46 -645827356227815.0,-1146491233835383.0,-1551601459626981.0,-1908764696253222.0,
47 -2234961611612173.0,-2538787015856429.0,-2825360342097020.0,
48 -3098043823061022.0,-3359196018589429.0,-3610552233837226.0,
49 }};
50
51 void
_arb_hypgeom_airy_zero(arb_t res,const fmpz_t n,int which,slong prec)52 _arb_hypgeom_airy_zero(arb_t res, const fmpz_t n, int which, slong prec)
53 {
54 slong asymp_accuracy, wp;
55
56 if (fmpz_cmp_ui(n, 10) <= 0)
57 {
58 if (fmpz_sgn(n) <= 0)
59 {
60 flint_printf("Airy zero only defined for index >= 1\n");
61 flint_abort();
62 }
63
64 /* The asymptotic expansions work well except when n == 1, so
65 use precomputed starting intervals (also for the first
66 few larger n as a small optimization). */
67 arf_set_d(arb_midref(res), ldexp(initial[which][fmpz_get_ui(n)-1], -48));
68 mag_set_d(arb_radref(res), ldexp(1.0, -48));
69 asymp_accuracy = 48;
70 }
71 else
72 {
73 arb_t z, u, u2, u4, s;
74 fmpz_t c;
75
76 arb_init(z);
77 arb_init(u);
78 arb_init(u2);
79 arb_init(u4);
80 arb_init(s);
81 fmpz_init(c);
82
83 if (which == AI || which == BI_PRIME)
84 asymp_accuracy = 13 + 10 * (fmpz_bits(n) - 1);
85 else
86 {
87 fmpz_sub_ui(c, n, 1);
88 asymp_accuracy = 13 + 10 * (fmpz_bits(c) - 1);
89 }
90
91 wp = asymp_accuracy + 8;
92 /* Reduce precision since we may not need to do any Newton steps. */
93 if (which == AI || which == BI)
94 wp = FLINT_MIN(wp, prec + 8);
95
96 arb_const_pi(z, wp);
97 fmpz_mul_2exp(c, n, 2);
98 if (which == AI || which == BI_PRIME)
99 fmpz_sub_ui(c, c, 1);
100 else
101 fmpz_sub_ui(c, c, 3);
102 fmpz_mul_ui(c, c, 3);
103 arb_mul_fmpz(z, z, c, wp);
104 arb_mul_2exp_si(z, z, -3);
105
106 arb_inv(u, z, wp);
107 arb_mul(u2, u, u, wp);
108 arb_mul(u4, u2, u2, wp);
109
110 if (which == AI || which == BI)
111 {
112 /* u^8 truncation gives lower bound */
113 arb_mul_si(s, u4, -108056875, wp);
114 arb_addmul_si(s, u2, 6478500, wp);
115 arb_add_si(s, s, -967680, wp);
116 arb_mul(s, s, u4, wp);
117 arb_addmul_si(s, u2, 725760, wp);
118 arb_div_ui(s, s, 6967296, wp);
119
120 /* u^10 term gives upper bound */
121 arb_mul(u4, u4, u4, 10);
122 arb_mul(u4, u4, u2, 10);
123 arb_mul_ui(u4, u4, 486, 10);
124 }
125 else
126 {
127 /* u^8 truncation gives upper bound */
128 arb_mul_si(s, u4, 18683371, wp);
129 arb_addmul_si(s, u2, -1087338, wp);
130 arb_add_si(s, s, 151200, wp);
131 arb_mul(s, s, u4, wp);
132 arb_addmul_si(s, u2, -181440, wp);
133 arb_div_ui(s, s, 1244160, wp);
134
135 /* u^10 term gives lower bound */
136 arb_mul(u4, u4, u4, 10);
137 arb_mul(u4, u4, u2, 10);
138 arb_mul_ui(u4, u4, 477, 10);
139 arb_neg(u4, u4);
140 }
141
142 arb_mul_2exp_si(u4, u4, -1);
143 arb_add(s, s, u4, wp);
144 arb_add_error(s, u4);
145
146 arb_add_ui(s, s, 1, wp);
147 arb_root_ui(z, z, 3, wp);
148 arb_mul(z, z, z, wp);
149 arb_mul(res, z, s, wp);
150
151 arb_neg(res, res);
152
153 arb_clear(z);
154 arb_clear(u);
155 arb_clear(u2);
156 arb_clear(u4);
157 arb_clear(s);
158 fmpz_clear(c);
159 }
160
161 /* Do interval Newton steps for refinement. Important: for the
162 primed zeros, we need to do at least one interval Newton step to
163 validate the initial (tentative) inclusion. */
164 if (asymp_accuracy < prec || (which == AI_PRIME || which == BI_PRIME))
165 {
166 arb_t f, fprime, root;
167 mag_t C, r;
168 slong * steps;
169 slong step, extraprec;
170
171 arb_init(f);
172 arb_init(fprime);
173 arb_init(root);
174 mag_init(C);
175 mag_init(r);
176 steps = flint_malloc(sizeof(slong) * FLINT_BITS);
177
178 extraprec = 0.25 * fmpz_bits(n) + 8;
179 wp = asymp_accuracy + extraprec;
180
181 /* C = |f''| or |f'''| on the initial interval given by res */
182 /* f''(x) = xf(x) */
183 /* f'''(x) = xf'(x) + f(x) */
184 if (which == AI || which == AI_PRIME)
185 arb_hypgeom_airy(f, fprime, NULL, NULL, res, wp);
186 else
187 arb_hypgeom_airy(NULL, NULL, f, fprime, res, wp);
188
189 if (which == AI || which == BI)
190 arb_mul(f, f, res, wp);
191 else
192 arb_addmul(f, fprime, res, wp);
193
194 arb_get_mag(C, f);
195
196 step = 0;
197 steps[step] = prec;
198
199 while (steps[step] / 2 > asymp_accuracy - extraprec)
200 {
201 steps[step + 1] = steps[step] / 2;
202 step++;
203 }
204
205 arb_set(root, res);
206
207 for ( ; step >= 0; step--)
208 {
209 wp = steps[step] + extraprec;
210 wp = FLINT_MAX(wp, arb_rel_accuracy_bits(root) + extraprec);
211
212 /* store radius, set root to the midpoint */
213 mag_set(r, arb_radref(root));
214 mag_zero(arb_radref(root));
215
216 if (which == AI || which == AI_PRIME)
217 arb_hypgeom_airy(f, fprime, NULL, NULL, root, wp);
218 else
219 arb_hypgeom_airy(NULL, NULL, f, fprime, root, wp);
220
221 /* f, f' = f', xf */
222 if (which == AI_PRIME || which == BI_PRIME)
223 {
224 arb_mul(f, f, root, wp);
225 arb_swap(f, fprime);
226 }
227
228 /* f'([m+/-r]) = f'(m) +/- f''([m +/- r]) * r */
229 mag_mul(r, C, r);
230 arb_add_error_mag(fprime, r);
231 arb_div(f, f, fprime, wp);
232 arb_sub(root, root, f, wp);
233
234 /* Verify inclusion so that C is still valid, and for the
235 primed zeros also to make sure that the initial
236 intervals really were correct. */
237 if (!arb_contains(res, root))
238 {
239 flint_printf("unexpected: no containment of Airy zero\n");
240 arb_indeterminate(root);
241 break;
242 }
243 }
244
245 arb_set(res, root);
246
247 arb_clear(f);
248 arb_clear(fprime);
249 arb_clear(root);
250 mag_clear(C);
251 mag_clear(r);
252 flint_free(steps);
253 }
254
255 arb_set_round(res, res, prec);
256 }
257
258 void
arb_hypgeom_airy_zero(arb_t ai,arb_t aip,arb_t bi,arb_t bip,const fmpz_t n,slong prec)259 arb_hypgeom_airy_zero(arb_t ai, arb_t aip, arb_t bi, arb_t bip, const fmpz_t n, slong prec)
260 {
261 if (ai != NULL)
262 _arb_hypgeom_airy_zero(ai, n, AI, prec);
263 if (aip != NULL)
264 _arb_hypgeom_airy_zero(aip, n, AI_PRIME, prec);
265 if (bi != NULL)
266 _arb_hypgeom_airy_zero(bi, n, BI, prec);
267 if (bip != NULL)
268 _arb_hypgeom_airy_zero(bip, n, BI_PRIME, prec);
269 }
270