mirror of
https://github.com/vale981/arb
synced 2025-03-05 09:21:38 -05:00
Computation of Airy function zeros
This commit is contained in:
fredrik
parent
492cd192bf
commit
d30df8c5a4
5 changed files with 452 additions and 0 deletions
|
|
@ -84,6 +84,8 @@ void arb_hypgeom_airy_series(arb_poly_t ai, arb_poly_t ai_prime,
|
|||
void _arb_hypgeom_airy_series(arb_ptr ai, arb_ptr ai_prime,
|
||||
arb_ptr bi, arb_ptr bi_prime, arb_srcptr z, slong zlen, slong len, slong prec);
|
||||
|
||||
void arb_hypgeom_airy_zero(arb_t ai, arb_t aip, arb_t bi, arb_t bip, const fmpz_t n, slong prec);
|
||||
|
||||
void arb_hypgeom_expint(arb_t res, const arb_t s, const arb_t z, slong prec);
|
||||
|
||||
void arb_hypgeom_gamma_lower(arb_t res, const arb_t s, const arb_t z, int regularized, slong prec);
|
||||
|
|
|
|||
269
arb_hypgeom/airy_zero.c
Normal file
269
arb_hypgeom/airy_zero.c
Normal file
|
|
@ -0,0 +1,269 @@
|
|||
/*
|
||||
Copyright (C) 2018 Fredrik Johansson
|
||||
|
||||
This file is part of Arb.
|
||||
|
||||
Arb is free software: you can redistribute it and/or modify it under
|
||||
the terms of the GNU Lesser General Public License (LGPL) as published
|
||||
by the Free Software Foundation; either version 2.1 of the License, or
|
||||
(at your option) any later version. See <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#include "arb_hypgeom.h"
|
||||
|
||||
/*
|
||||
https://dlmf.nist.gov/9.9
|
||||
|
||||
a_k ~ -T(3/8 pi (4k-1))
|
||||
a'_k ~ -U(3/8 pi (4k-3))
|
||||
b_k ~ -T(3/8 pi (4k-3))
|
||||
b'_k ~ -U(3/8 pi (4k-1))
|
||||
|
||||
For a_k and b_k, the u^8 and u^10 truncations are known to give lower
|
||||
bounds. [G. Pittaluga and L. Sacripante (1991) Inequalities for the
|
||||
zeros of the Airy functions. SIAM J. Math. Anal. 22 (1), pp. 260–267.]
|
||||
|
||||
We don't have proofs for a'_k and b'_k. However, in that case, we can just
|
||||
do a single interval Newton step to verify that we have isolated a
|
||||
zero (the enclosure must be for the correct zero due to sandwiching).
|
||||
*/
|
||||
|
||||
#define AI 0
|
||||
#define BI 1
|
||||
#define AI_PRIME 2
|
||||
#define BI_PRIME 3
|
||||
|
||||
static const double initial[4][10] = {{
|
||||
-658118728906175.0,-1150655474581104.0,-1553899449042978.0,-1910288501594969.0,
|
||||
-2236074816421182.0,-2539650438812533.0,-2826057838960988.0,
|
||||
-3098624122012011.0,-3359689702679955.0,-3610979637739094.0},{
|
||||
-330370902027041.0,-920730911234245.0,-1359731821477101.0,-1736658984124319.0,
|
||||
-2076373934490092.0,-2390271103799312.0,-2684763040788193.0,
|
||||
-2963907159065113.0,-3230475233555475.0,-3486466475611047.0},{
|
||||
-286764727967452.0,-914286338795679.0,-1356737313209586.0,-1734816794389239.0,
|
||||
-2075083421171399.0,-2389296605766914.0,-2683990299959380.0,
|
||||
-2963272965051282.0,-3229941298662311.0,-3486008018531685.0},{
|
||||
-645827356227815.0,-1146491233835383.0,-1551601459626981.0,-1908764696253222.0,
|
||||
-2234961611612173.0,-2538787015856429.0,-2825360342097020.0,
|
||||
-3098043823061022.0,-3359196018589429.0,-3610552233837226.0,
|
||||
}};
|
||||
|
||||
void
|
||||
_arb_hypgeom_airy_zero(arb_t res, const fmpz_t n, int which, slong prec)
|
||||
{
|
||||
slong asymp_accuracy, wp;
|
||||
|
||||
if (fmpz_cmp_ui(n, 10) <= 0)
|
||||
{
|
||||
if (fmpz_sgn(n) <= 0)
|
||||
{
|
||||
flint_printf("Airy zero only defined for index >= 1\n");
|
||||
flint_abort();
|
||||
}
|
||||
|
||||
/* The asymptotic expansions work well except when n == 1, so
|
||||
use precomputed starting intervals (also for the first
|
||||
few larger n as a small optimization). */
|
||||
arf_set_d(arb_midref(res), ldexp(initial[which][fmpz_get_ui(n)-1], -48));
|
||||
mag_set_d(arb_radref(res), ldexp(1.0, -48));
|
||||
asymp_accuracy = 48;
|
||||
}
|
||||
else
|
||||
{
|
||||
arb_t z, u, u2, u4, s;
|
||||
fmpz_t c;
|
||||
|
||||
arb_init(z);
|
||||
arb_init(u);
|
||||
arb_init(u2);
|
||||
arb_init(u4);
|
||||
arb_init(s);
|
||||
fmpz_init(c);
|
||||
|
||||
if (which == AI || which == BI_PRIME)
|
||||
asymp_accuracy = 13 + 10 * (fmpz_bits(n) - 1);
|
||||
else
|
||||
{
|
||||
fmpz_sub_ui(c, n, 1);
|
||||
asymp_accuracy = 13 + 10 * (fmpz_bits(c) - 1);
|
||||
}
|
||||
|
||||
wp = asymp_accuracy + 8;
|
||||
/* Reduce precision since we may not need to do any Newton steps. */
|
||||
if (which == AI || which == BI)
|
||||
wp = FLINT_MIN(wp, prec + 8);
|
||||
|
||||
arb_const_pi(z, wp);
|
||||
fmpz_mul_2exp(c, n, 2);
|
||||
if (which == AI || which == BI_PRIME)
|
||||
fmpz_sub_ui(c, c, 1);
|
||||
else
|
||||
fmpz_sub_ui(c, c, 3);
|
||||
fmpz_mul_ui(c, c, 3);
|
||||
arb_mul_fmpz(z, z, c, wp);
|
||||
arb_mul_2exp_si(z, z, -3);
|
||||
|
||||
arb_inv(u, z, wp);
|
||||
arb_mul(u2, u, u, wp);
|
||||
arb_mul(u4, u2, u2, wp);
|
||||
|
||||
if (which == AI || which == BI)
|
||||
{
|
||||
/* u^8 truncation gives lower bound */
|
||||
arb_mul_si(s, u4, -108056875, wp);
|
||||
arb_addmul_si(s, u2, 6478500, wp);
|
||||
arb_add_si(s, s, -967680, wp);
|
||||
arb_mul(s, s, u4, wp);
|
||||
arb_addmul_si(s, u2, 725760, wp);
|
||||
arb_div_ui(s, s, 6967296, wp);
|
||||
|
||||
/* u^10 term gives upper bound */
|
||||
arb_mul(u4, u4, u4, 10);
|
||||
arb_mul(u4, u4, u2, 10);
|
||||
arb_mul_ui(u4, u4, 486, 10);
|
||||
}
|
||||
else
|
||||
{
|
||||
/* u^8 truncation gives upper bound */
|
||||
arb_mul_si(s, u4, 18683371, wp);
|
||||
arb_addmul_si(s, u2, -1087338, wp);
|
||||
arb_add_si(s, s, 151200, wp);
|
||||
arb_mul(s, s, u4, wp);
|
||||
arb_addmul_si(s, u2, -181440, wp);
|
||||
arb_div_ui(s, s, 1244160, wp);
|
||||
|
||||
/* u^10 term gives lower bound */
|
||||
arb_mul(u4, u4, u4, 10);
|
||||
arb_mul(u4, u4, u2, 10);
|
||||
arb_mul_ui(u4, u4, 477, 10);
|
||||
arb_neg(u4, u4);
|
||||
}
|
||||
|
||||
arb_mul_2exp_si(u4, u4, -1);
|
||||
arb_add(s, s, u4, wp);
|
||||
arb_add_error(s, u4);
|
||||
|
||||
arb_add_ui(s, s, 1, wp);
|
||||
arb_root_ui(z, z, 3, wp);
|
||||
arb_mul(z, z, z, wp);
|
||||
arb_mul(res, z, s, wp);
|
||||
|
||||
arb_neg(res, res);
|
||||
|
||||
arb_clear(z);
|
||||
arb_clear(u);
|
||||
arb_clear(u2);
|
||||
arb_clear(u4);
|
||||
arb_clear(s);
|
||||
fmpz_clear(c);
|
||||
}
|
||||
|
||||
/* Do interval Newton steps for refinement. Important: for the
|
||||
primed zeros, we need to do at least one interval Newton step to
|
||||
validate the initial (tentative) inclusion. */
|
||||
if (asymp_accuracy < prec || (which == AI_PRIME || which == BI_PRIME))
|
||||
{
|
||||
arb_t f, fprime, root;
|
||||
mag_t C, r;
|
||||
slong * steps;
|
||||
slong step, extraprec;
|
||||
|
||||
arb_init(f);
|
||||
arb_init(fprime);
|
||||
arb_init(root);
|
||||
mag_init(C);
|
||||
mag_init(r);
|
||||
steps = flint_malloc(sizeof(slong) * FLINT_BITS);
|
||||
|
||||
extraprec = 0.25 * fmpz_bits(n) + 8;
|
||||
wp = asymp_accuracy + extraprec;
|
||||
|
||||
/* C = |f''| or |f'''| on the initial interval given by res */
|
||||
/* f''(x) = xf(x) */
|
||||
/* f'''(x) = xf'(x) + f(x) */
|
||||
if (which == AI || which == AI_PRIME)
|
||||
arb_hypgeom_airy(f, fprime, NULL, NULL, res, wp);
|
||||
else
|
||||
arb_hypgeom_airy(NULL, NULL, f, fprime, res, wp);
|
||||
|
||||
if (which == AI || which == BI)
|
||||
arb_mul(f, f, res, wp);
|
||||
else
|
||||
arb_addmul(f, fprime, res, wp);
|
||||
|
||||
arb_get_mag(C, f);
|
||||
|
||||
step = 0;
|
||||
steps[step] = prec;
|
||||
|
||||
while (steps[step] / 2 > asymp_accuracy - extraprec)
|
||||
{
|
||||
steps[step + 1] = steps[step] / 2;
|
||||
step++;
|
||||
}
|
||||
|
||||
arb_set(root, res);
|
||||
|
||||
for ( ; step >= 0; step--)
|
||||
{
|
||||
wp = steps[step] + extraprec;
|
||||
wp = FLINT_MAX(wp, arb_rel_accuracy_bits(root) + extraprec);
|
||||
|
||||
/* store radius, set root to the midpoint */
|
||||
mag_set(r, arb_radref(root));
|
||||
mag_zero(arb_radref(root));
|
||||
|
||||
if (which == AI || which == AI_PRIME)
|
||||
arb_hypgeom_airy(f, fprime, NULL, NULL, root, wp);
|
||||
else
|
||||
arb_hypgeom_airy(NULL, NULL, f, fprime, root, wp);
|
||||
|
||||
/* f, f' = f', xf */
|
||||
if (which == AI_PRIME || which == BI_PRIME)
|
||||
{
|
||||
arb_mul(f, f, root, wp);
|
||||
arb_swap(f, fprime);
|
||||
}
|
||||
|
||||
/* f'([m+/-r]) = f'(m) +/- f''([m +/- r]) * r */
|
||||
mag_mul(r, C, r);
|
||||
arb_add_error_mag(fprime, r);
|
||||
arb_div(f, f, fprime, wp);
|
||||
arb_sub(root, root, f, wp);
|
||||
|
||||
/* Verify inclusion so that C is still valid, and for the
|
||||
primed zeros also to make sure that the initial
|
||||
intervals really were correct. */
|
||||
if (!arb_contains(res, root))
|
||||
{
|
||||
flint_printf("unexpected: no containment of Airy zero\n");
|
||||
arb_indeterminate(root);
|
||||
break;
|
||||
}
|
||||
}
|
||||
|
||||
arb_set(res, root);
|
||||
|
||||
arb_clear(f);
|
||||
arb_clear(fprime);
|
||||
arb_clear(root);
|
||||
mag_clear(C);
|
||||
mag_clear(r);
|
||||
flint_free(steps);
|
||||
}
|
||||
|
||||
arb_set_round(res, res, prec);
|
||||
}
|
||||
|
||||
void
|
||||
arb_hypgeom_airy_zero(arb_t ai, arb_t aip, arb_t bi, arb_t bip, const fmpz_t n, slong prec)
|
||||
{
|
||||
if (ai != NULL)
|
||||
_arb_hypgeom_airy_zero(ai, n, AI, prec);
|
||||
if (aip != NULL)
|
||||
_arb_hypgeom_airy_zero(aip, n, AI_PRIME, prec);
|
||||
if (bi != NULL)
|
||||
_arb_hypgeom_airy_zero(bi, n, BI, prec);
|
||||
if (bip != NULL)
|
||||
_arb_hypgeom_airy_zero(bip, n, BI_PRIME, prec);
|
||||
}
|
||||
168
arb_hypgeom/test/t-airy_zero.c
Normal file
168
arb_hypgeom/test/t-airy_zero.c
Normal file
|
|
@ -0,0 +1,168 @@
|
|||
/*
|
||||
Copyright (C) 2018 Fredrik Johansson
|
||||
|
||||
This file is part of Arb.
|
||||
|
||||
Arb is free software: you can redistribute it and/or modify it under
|
||||
the terms of the GNU Lesser General Public License (LGPL) as published
|
||||
by the Free Software Foundation; either version 2.1 of the License, or
|
||||
(at your option) any later version. See <http://www.gnu.org/licenses/>.
|
||||
*/
|
||||
|
||||
#include "arb_hypgeom.h"
|
||||
|
||||
int main()
|
||||
{
|
||||
slong iter;
|
||||
flint_rand_t state;
|
||||
|
||||
flint_printf("airy_zero....");
|
||||
fflush(stdout);
|
||||
|
||||
flint_randinit(state);
|
||||
|
||||
/* interlacing test */
|
||||
{
|
||||
arb_t a, ap, b, bp, ap1;
|
||||
fmpz_t n;
|
||||
|
||||
arb_init(a);
|
||||
arb_init(ap);
|
||||
arb_init(b);
|
||||
arb_init(bp);
|
||||
arb_init(ap1);
|
||||
fmpz_init(n);
|
||||
|
||||
for (fmpz_one(n); fmpz_cmp_ui(n, 200) <= 0; fmpz_add_ui(n, n, 1))
|
||||
{
|
||||
arb_hypgeom_airy_zero(a, ap, b, bp, n, 53);
|
||||
fmpz_add_ui(n, n, 1);
|
||||
arb_hypgeom_airy_zero(NULL, ap1, NULL, NULL, n, 53);
|
||||
fmpz_sub_ui(n, n, 1);
|
||||
|
||||
if (!arb_gt(ap, b) || !arb_gt(b, bp) || !arb_gt(bp, a) || !arb_gt(a, ap1))
|
||||
{
|
||||
flint_printf("FAIL: interlacing\n\n");
|
||||
flint_printf("n = "); fmpz_print(n); flint_printf("\n\n");
|
||||
flint_printf("a = "); arb_printn(a, 100, 0); flint_printf("\n\n");
|
||||
flint_printf("ap = "); arb_printn(ap, 100, 0); flint_printf("\n\n");
|
||||
flint_printf("b = "); arb_printn(b, 100, 0); flint_printf("\n\n");
|
||||
flint_printf("bp = "); arb_printn(bp, 100, 0); flint_printf("\n\n");
|
||||
flint_printf("ap1 = "); arb_printn(ap1, 100, 0); flint_printf("\n\n");
|
||||
flint_abort();
|
||||
}
|
||||
}
|
||||
|
||||
for (fmpz_set_ui(n, 1000000); fmpz_cmp_ui(n, 1000000 + 200) <= 0; fmpz_add_ui(n, n, 1))
|
||||
{
|
||||
arb_hypgeom_airy_zero(a, ap, b, bp, n, 53);
|
||||
fmpz_add_ui(n, n, 1);
|
||||
arb_hypgeom_airy_zero(NULL, ap1, NULL, NULL, n, 53);
|
||||
fmpz_sub_ui(n, n, 1);
|
||||
|
||||
if (!arb_gt(ap, b) || !arb_gt(b, bp) || !arb_gt(bp, a) || !arb_gt(a, ap1))
|
||||
{
|
||||
flint_printf("FAIL: interlacing\n\n");
|
||||
flint_printf("n = "); fmpz_print(n); flint_printf("\n\n");
|
||||
flint_printf("a = "); arb_printn(a, 100, 0); flint_printf("\n\n");
|
||||
flint_printf("ap = "); arb_printn(ap, 100, 0); flint_printf("\n\n");
|
||||
flint_printf("b = "); arb_printn(b, 100, 0); flint_printf("\n\n");
|
||||
flint_printf("bp = "); arb_printn(bp, 100, 0); flint_printf("\n\n");
|
||||
flint_printf("ap1 = "); arb_printn(ap1, 100, 0); flint_printf("\n\n");
|
||||
flint_abort();
|
||||
}
|
||||
}
|
||||
|
||||
arb_clear(a);
|
||||
arb_clear(ap);
|
||||
arb_clear(b);
|
||||
arb_clear(bp);
|
||||
arb_clear(ap1);
|
||||
fmpz_clear(n);
|
||||
}
|
||||
|
||||
/* self-consistency and accuracy test */
|
||||
for (iter = 0; iter < 2000 * arb_test_multiplier(); iter++)
|
||||
{
|
||||
arb_t x1, x2, v1, v2;
|
||||
fmpz_t n;
|
||||
slong prec1, prec2;
|
||||
int which;
|
||||
|
||||
arb_init(x1);
|
||||
arb_init(x2);
|
||||
arb_init(v1);
|
||||
arb_init(v2);
|
||||
fmpz_init(n);
|
||||
|
||||
fmpz_randtest(n, state, 200);
|
||||
fmpz_abs(n, n);
|
||||
fmpz_add_ui(n, n, 1);
|
||||
prec1 = 2 + n_randtest(state) % 1000;
|
||||
prec2 = 2 + n_randtest(state) % 1000;
|
||||
which = n_randint(state, 4);
|
||||
|
||||
if (which == 0)
|
||||
{
|
||||
arb_hypgeom_airy_zero(x1, NULL, NULL, NULL, n, prec1);
|
||||
arb_hypgeom_airy_zero(x2, NULL, NULL, NULL, n, prec2);
|
||||
arb_hypgeom_airy(v1, NULL, NULL, NULL, x1, prec1 + 30);
|
||||
arb_hypgeom_airy(v2, NULL, NULL, NULL, x2, prec2 + 30);
|
||||
}
|
||||
else if (which == 1)
|
||||
{
|
||||
arb_hypgeom_airy_zero(NULL, x1, NULL, NULL, n, prec1);
|
||||
arb_hypgeom_airy_zero(NULL, x2, NULL, NULL, n, prec2);
|
||||
arb_hypgeom_airy(NULL, v1, NULL, NULL, x1, prec1 + 30);
|
||||
arb_hypgeom_airy(NULL, v2, NULL, NULL, x2, prec2 + 30);
|
||||
}
|
||||
else if (which == 2)
|
||||
{
|
||||
arb_hypgeom_airy_zero(NULL, NULL, x1, NULL, n, prec1);
|
||||
arb_hypgeom_airy_zero(NULL, NULL, x2, NULL, n, prec2);
|
||||
arb_hypgeom_airy(NULL, NULL, v1, NULL, x1, prec1 + 30);
|
||||
arb_hypgeom_airy(NULL, NULL, v2, NULL, x2, prec2 + 30);
|
||||
}
|
||||
else
|
||||
{
|
||||
arb_hypgeom_airy_zero(NULL, NULL, NULL, x1, n, prec1);
|
||||
arb_hypgeom_airy_zero(NULL, NULL, NULL, x2, n, prec2);
|
||||
arb_hypgeom_airy(NULL, NULL, NULL, v1, x1, prec1 + 30);
|
||||
arb_hypgeom_airy(NULL, NULL, NULL, v2, x2, prec2 + 30);
|
||||
}
|
||||
|
||||
if (!arb_overlaps(x1, x2) || !arb_contains_zero(v1) || !arb_contains_zero(v2))
|
||||
{
|
||||
flint_printf("FAIL: overlap\n\n");
|
||||
flint_printf("which = %d, n = ", which); fmpz_print(n);
|
||||
flint_printf(" prec1 = %wd prec2 = %wd\n\n", prec1, prec2);
|
||||
flint_printf("x1 = "); arb_printn(x1, 100, 0); flint_printf("\n\n");
|
||||
flint_printf("x2 = "); arb_printn(x2, 100, 0); flint_printf("\n\n");
|
||||
flint_printf("v1 = "); arb_printn(v1, 100, 0); flint_printf("\n\n");
|
||||
flint_printf("v2 = "); arb_printn(v2, 100, 0); flint_printf("\n\n");
|
||||
flint_abort();
|
||||
}
|
||||
|
||||
if (arb_rel_accuracy_bits(x1) < prec1 - 3 || arb_rel_accuracy_bits(x2) < prec2 - 3)
|
||||
{
|
||||
flint_printf("FAIL: accuracy\n\n");
|
||||
flint_printf("which = %d, n = ", which); fmpz_print(n);
|
||||
flint_printf(" prec1 = %wd prec2 = %wd\n\n", prec1, prec2);
|
||||
flint_printf("acc(x1) = %wd, acc(x2) = %wd\n\n", arb_rel_accuracy_bits(x1), arb_rel_accuracy_bits(x2));
|
||||
flint_printf("x1 = "); arb_printn(x1, 100, 0); flint_printf("\n\n");
|
||||
flint_printf("x2 = "); arb_printn(x2, 100, 0); flint_printf("\n\n");
|
||||
flint_abort();
|
||||
}
|
||||
|
||||
arb_clear(x1);
|
||||
arb_clear(x2);
|
||||
arb_clear(v1);
|
||||
arb_clear(v2);
|
||||
fmpz_clear(n);
|
||||
}
|
||||
|
||||
flint_randclear(state);
|
||||
flint_cleanup();
|
||||
flint_printf("PASS\n");
|
||||
return EXIT_SUCCESS;
|
||||
}
|
||||
|
|
@ -322,6 +322,17 @@ Airy functions
|
|||
truncated to length *len*. As with the other Airy methods, any of the
|
||||
outputs can be *NULL*.
|
||||
|
||||
.. function:: void arb_hypgeom_airy_zero(arb_t a, arb_t a_prime, arb_t b, arb_t b_prime, const fmpz_t n, slong prec)
|
||||
|
||||
Computes the *n*-th real zero `a_n`, `a'_n`, `b_n`, or `b'_n`
|
||||
for the respective Airy function or Airy function derivative.
|
||||
Any combination of the four output variables can be *NULL*.
|
||||
The zeros are indexed by increasing magnitude, starting with
|
||||
`n = 1` to follow the convention in the literature.
|
||||
An index *n* that is not positive is invalid input.
|
||||
The implementation uses asymptotic expansions for the zeros
|
||||
[PS1991]_ together with the interval Newton method for refinement.
|
||||
|
||||
Orthogonal polynomials and functions
|
||||
-------------------------------------------------------------------------------
|
||||
|
||||
|
|
|
|||
|
|
@ -236,6 +236,8 @@ Bibliography
|
|||
|
||||
.. [PS1973] \M. S. Paterson and L. J. Stockmeyer, "On the number of nonscalar multiplications necessary to evaluate polynomials", SIAM J. Comput (1973)
|
||||
|
||||
.. [PS1991] \G. Pittaluga and L. Sacripante, "Inequalities for the zeros of the Airy functions", SIAM J. Math. Anal. 22:1 (1991), 260-267.
|
||||
|
||||
.. [Rum2010] \S. M. Rump, "Verification methods: Rigorous results using floating-point arithmetic", Acta Numerica 19 (2010), 287-449.
|
||||
|
||||
.. [Smi2001] \D. M. Smith, "Algorithm: Fortran 90 Software for Floating-Point Multiple Precision Arithmetic, Gamma and Related Functions", Transactions on Mathematical Software 27 (2001) 377-387, http://myweb.lmu.edu/dmsmith/toms2001.pdf
|
||||
|
|
|
|||
Loading…
Add table
Reference in a new issue