diff --git a/arb_hypgeom.h b/arb_hypgeom.h
index cadc95bd..1dd93f75 100644
--- a/arb_hypgeom.h
+++ b/arb_hypgeom.h
@@ -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);
diff --git a/arb_hypgeom/airy_zero.c b/arb_hypgeom/airy_zero.c
new file mode 100644
index 00000000..73eaf948
--- /dev/null
+++ b/arb_hypgeom/airy_zero.c
@@ -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 .
+*/
+
+#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);
+}
diff --git a/arb_hypgeom/test/t-airy_zero.c b/arb_hypgeom/test/t-airy_zero.c
new file mode 100644
index 00000000..f1c84401
--- /dev/null
+++ b/arb_hypgeom/test/t-airy_zero.c
@@ -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 .
+*/
+
+#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;
+}
diff --git a/doc/source/arb_hypgeom.rst b/doc/source/arb_hypgeom.rst
index b7faf0c7..77713834 100644
--- a/doc/source/arb_hypgeom.rst
+++ b/doc/source/arb_hypgeom.rst
@@ -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
-------------------------------------------------------------------------------
diff --git a/doc/source/credits.rst b/doc/source/credits.rst
index d151900b..8d4430b8 100644
--- a/doc/source/credits.rst
+++ b/doc/source/credits.rst
@@ -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