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