attempt to improve bounds for Airy functions

acb_hypgeom/airy.c

@@ -102,8 +102,8 @@ estimate_airy(double x, double y, int ai)
 
 /* error propagation based on derivatives */
 void
-acb_hypgeom_airy_direct_prop(acb_t ai, acb_t aip, acb_t bi, acb_t bip,
+acb_hypgeom_airy_prop(acb_t ai, acb_t aip, acb_t bi, acb_t bip,
-    const acb_t z, slong n, slong prec)
+    const acb_t z, slong n, int algo, slong prec)
 {
     mag_t aib, aipb, bib, bipb, zb, rad;
     acb_t zz;
@@ -124,7 +124,10 @@ acb_hypgeom_airy_direct_prop(acb_t ai, acb_t aip, acb_t bi, acb_t bip,
     acb_get_mag(zb, z);
 
     acb_hypgeom_airy_bound(aib, aipb, bib, bipb, z);
+    if (algo == 0)
         acb_hypgeom_airy_direct(ai, aip, bi, bip, zz, n, prec);
+    else
+        acb_hypgeom_airy_asymp(ai, aip, bi, bip, zz, n, prec);
 
     if (ai != NULL)
     {
@@ -173,6 +176,24 @@ acb_hypgeom_airy_direct_prop(acb_t ai, acb_t aip, acb_t bi, acb_t bip,
     acb_clear(zz);
 }
 
+void
+acb_hypgeom_airy_direct_prop(acb_t ai, acb_t aip, acb_t bi, acb_t bip,
+    const acb_t z, slong n, slong prec)
+{
+    acb_hypgeom_airy_prop(ai, aip, bi, bip, z, n, 0, prec);
+}
+
+void
+acb_hypgeom_airy_asymp2(acb_t ai, acb_t aip, acb_t bi, acb_t bip,
+    const acb_t z, slong n, slong prec)
+{
+    /* avoid singularity in asymptotic expansion near 0 */
+    if (acb_rel_accuracy_bits(z) > 3)
+        acb_hypgeom_airy_asymp(ai, aip, bi, bip, z, n, prec);
+    else
+        acb_hypgeom_airy_prop(ai, aip, bi, bip, z, n, 1, prec);
+}
+
 void
 acb_hypgeom_airy(acb_t ai, acb_t aip, acb_t bi, acb_t bip, const acb_t z, slong prec)
 {
@@ -210,7 +231,10 @@ acb_hypgeom_airy(acb_t ai, acb_t aip, acb_t bi, acb_t bip, const acb_t z, slong
             n = wp / (-zmag) + 1;
         }
 
+        if (acb_is_exact(z))
             acb_hypgeom_airy_direct(ai, aip, bi, bip, z, n, wp);
+        else
+            acb_hypgeom_airy_direct_prop(ai, aip, bi, bip, z, n, wp);
     }  /* huge input -- use asymptotics and pick n without overflowing */
     else if ((arf_cmpabs_2exp_si(re, 64) > 0 || arf_cmpabs_2exp_si(im, 64) > 0))
     {
@@ -227,7 +251,7 @@ acb_hypgeom_airy(acb_t ai, acb_t aip, acb_t bi, acb_t bip, const acb_t z, slong
             n = FLINT_MAX(n, 1);
         }
 
-        acb_hypgeom_airy_asymp(ai, aip, bi, bip, z, n, wp);
+        acb_hypgeom_airy_asymp2(ai, aip, bi, bip, z, n, wp);
     }
     else /* moderate input */
     {
@@ -239,14 +263,18 @@ acb_hypgeom_airy(acb_t ai, acb_t aip, acb_t bi, acb_t bip, const acb_t z, slong
 
         if (zmag >= 4.0 && (n = asymp_pick_terms(wp, log(zmag))) != -1)
         {
-            acb_hypgeom_airy_asymp(ai, aip, bi, bip, z, n, wp);
+            acb_hypgeom_airy_asymp2(ai, aip, bi, bip, z, n, wp);
         }
         else if (zmag <= 1.5)
         {
             t = 3 * (wp * LOG2) / (2 * z15 * EXP1);
             t = (wp * LOG2) / (2 * d_lambertw(t));
             n = FLINT_MAX(t + 1, 2);
+
+            if (acb_is_exact(z))
                 acb_hypgeom_airy_direct(ai, aip, bi, bip, z, n, wp);
+            else
+                acb_hypgeom_airy_direct_prop(ai, aip, bi, bip, z, n, wp);
         }
         else
         {

acb_hypgeom/airy_bound.c

@@ -24,65 +24,29 @@ arb_bound_exp_neg(mag_t b, const arb_t x)
     arb_clear(t);
 }
 
-/* todo -- should be lt in asymp code? */
+/* Implements DLMF 9.7.17. We assume |zeta| >= 1/2 and |arg(z)| <= 2 pi/3 here,
-static int
+   ignoring the smaller points which must be dealt with separately. */
-arg_le_2pi3(const acb_t z, const acb_t zeta)
-{
-    if (arb_is_nonnegative(acb_realref(z)))
-        return 1;
-
-    if (arb_is_positive(acb_imagref(z)) &&
-        arb_is_nonnegative(acb_imagref(zeta)))
-        return 1;
-
-    if (arb_is_negative(acb_imagref(z)) &&
-        arb_is_nonpositive(acb_imagref(zeta)))
-        return 1;
-
-    return 0;
-}
-
 void
 acb_hypgeom_airy_bound_9_7_17(mag_t bound, const acb_t z, const acb_t zeta)
 {
-    mag_t D, t, u, v, zeta_lower;
+    mag_t D, t, u, v, zeta_lower, half;
 
     mag_init(D);
     mag_init(t);
     mag_init(u);
     mag_init(v);
     mag_init(zeta_lower);
+    mag_init(half);
+
+    mag_one(half);
+    mag_mul_2exp_si(half, half, -1);
 
     acb_get_mag_lower(zeta_lower, zeta);
+    mag_max(zeta_lower, zeta_lower, half);  /* by assumption */
 
     /* 2 chi(1) exp(7 pi / (72 |zeta|)) * c_1 */
-    /* simplified bound */
+    /* simplified bound assuming |zeta| >= 1/2 */
-    if (mag_cmp_2exp_si(zeta_lower, -1) >= 0)
     mag_one(D);
-    else
-        mag_inf(D);
-
-    if (!arg_le_2pi3(z, zeta))
-    {
-        arb_get_mag_lower(u, acb_realref(zeta));
-        arb_get_mag(v, acb_imagref(zeta));
-
-        /* exp(7 pi / (36 u)) < exp(5/(8 u)) */
-        mag_set_ui_2exp_si(t, 5, -3);
-        mag_div(t, t, u);
-        mag_exp(t, t);
-
-        /* |1/cos(arg(zeta))| = sqrt(1+(v/u)^2) */
-        mag_div(v, v, u);
-        mag_mul(v, v, v);
-        mag_one(u);
-        mag_add(v, v, u);
-        mag_sqrt(v, v);
-        mag_mul(t, t, v);
-        /* c_1 * 4 chi(1) < 0.62 < 1 -- do nothing */
-
-        mag_max(D, D, t);
-    }
 
     /* exp(-zeta) / (2 sqrt(pi)) * (1 + D / |zeta|) */
 
@@ -106,38 +70,60 @@ acb_hypgeom_airy_bound_9_7_17(mag_t bound, const acb_t z, const acb_t zeta)
 }
 
 void
-acb_hypgeom_airy_bound(mag_t ai, mag_t aip, mag_t bi, mag_t bip, const acb_t z)
+acb_hypgeom_airy_bound_arg_le_2pi3(mag_t A, mag_t B, const acb_t z, slong wp)
 {
-    acb_t zeta;
+    acb_t zeta, z1;
-    acb_t z1, z2;
-    arf_srcptr zre, zim;
-    mag_t A, B, D, zlo, zhi;
-    slong wp;
-
-    if (acb_contains_zero(z))
-    {
-        if (ai != NULL) mag_inf(ai);
-        if (aip != NULL) mag_inf(aip);
-        if (bi != NULL) mag_inf(bi);
-        if (bip != NULL) mag_inf(bip);
-        return;
-    }
 
     acb_init(zeta);
-    mag_init(A);
+    acb_init(z1);
-    mag_init(B);
-    mag_init(D);
-    mag_init(zlo);
-    mag_init(zhi);
 
-    wp = MAG_BITS * 2;
+    acb_set_round(zeta, z, wp);
-    zre = arb_midref(acb_realref(z));
+    acb_sqrt(zeta, zeta, wp);
-    zim = arb_midref(acb_imagref(z));
+    acb_cube(zeta, zeta, wp);
+    acb_mul_2exp_si(zeta, zeta, 1);
+    acb_div_ui(zeta, zeta, 3, wp);
 
-    /* near the negative half line, use 
+    acb_hypgeom_airy_bound_9_7_17(A, z, zeta);
-       |Ai(-z)|, |Bi(-z)| <= |Ai(z exp(pi i/3))| + |Ai(z exp(-pi i/3))| */
+
-    if (arf_sgn(zre) < 0 && arf_cmpabs(zre, zim) > 0)
+    /* Use Bi(z) = w_1 Ai(z) + 2 w_2 Ai(z exp(+/- 2pi i / 3)),
+       where w_1, w_2 are roots of unity */
+    if (B != NULL)
     {
+        arb_sqrt_ui(acb_imagref(z1), 3, wp);
+        arb_set_si(acb_realref(z1), -1);
+        acb_mul_2exp_si(z1, z1, -1);
+
+        /* multiply by exp(-2 pi i / 3) in upper half plane
+           and by exp(2 pi i / 3) in lower half plane, to stay close
+           to positive reals */
+
+        /* in principle, we should be computing the union of both cases,
+           but since Bi is conjugate symmetric with the maximum value
+           attained on the positive real line, it's sufficient to consider
+           the case centered on the midpoint */
+        if (arf_sgn(arb_midref(acb_imagref(z))) >= 0)
+            acb_conj(z1, z1);
+
+        acb_mul(z1, z1, z, wp);
+        acb_neg(zeta, zeta);   /* same effect regardless of exp(+/-2 pi i/3) */
+
+        acb_hypgeom_airy_bound_9_7_17(B, z1, zeta);
+        mag_mul_2exp_si(B, B, 1);
+        mag_add(B, B, A);
+    }
+
+    acb_clear(zeta);
+    acb_clear(z1);
+}
+
+/* near the negative half line, use
+   |Ai(-z)|, |Bi(-z)| <= |Ai(z exp(pi i/3))| + |Ai(z exp(-pi i/3))| */
+void
+acb_hypgeom_airy_bound_arg_ge_2pi3(mag_t A, mag_t B, const acb_t z, slong wp)
+{
+    acb_t zeta, z1, z2;
+
+    acb_init(zeta);
     acb_init(z1);
     acb_init(z2);
 
@@ -164,56 +150,211 @@ acb_hypgeom_airy_bound(mag_t ai, mag_t aip, mag_t bi, mag_t bip, const acb_t z)
     }
     else
     {
+        mag_t D;
+        mag_init(D);
         acb_mul(z2, z2, zeta, wp);
         acb_neg(zeta, zeta);
         acb_hypgeom_airy_bound_9_7_17(D, z2, zeta);
         mag_add(A, A, D);
+        mag_clear(D);
     }
 
+    if (B != NULL)
         mag_set(B, A);
 
+    acb_clear(zeta);
     acb_clear(z1);
     acb_clear(z2);
+}
+
+static int
+arg_lt_2pi3_fast(const acb_t z)
+{
+    arf_t t;
+    mag_t x, y, s;
+    int res;
+
+    if (arb_is_zero(acb_imagref(z)) && arb_is_nonnegative(acb_realref(z)))
+        return 1;
+
+    arf_init(t);
+    mag_init(x);
+    mag_init(y);
+    mag_init(s);
+
+    arf_set_mag(t, arb_radref(acb_realref(z)));
+    arf_sub(t, arb_midref(acb_realref(z)), t, MAG_BITS, ARF_RND_FLOOR);
+
+    if (arf_sgn(t) >= 0)
+    {
+        res = 1;
     }
     else
     {
-        acb_set_round(zeta, z, wp);
+        arf_get_mag(x, t);
-        acb_sqrt(zeta, zeta, wp);
+        arb_get_mag_lower(y, acb_imagref(z));
-        acb_cube(zeta, zeta, wp);
+        mag_set_ui(s, 3);
-        acb_mul_2exp_si(zeta, zeta, 1);
+        mag_sqrt(s, s);
-        acb_div_ui(zeta, zeta, 3, wp);
+        mag_mul(s, s, x);
+        res = mag_cmp(s, y) <= 0;
+    }
 
-        acb_hypgeom_airy_bound_9_7_17(A, z, zeta);
+    arf_clear(t);
+    mag_clear(x);
+    mag_clear(y);
+    mag_clear(s);
 
-        /* Use Bi(z) = w_1 Ai(z) + 2 w_2 Ai(z exp(+/- 2pi i / 3)),
+    return res;
-           where w_1, w_2 are roots of unity */
+}
-        if (bi != NULL || bip != NULL)
+
+static int
+arg_gt_2pi3_fast(const acb_t z)
+{
+    arf_t t;
+    mag_t x, y, s;
+    int res;
+
+    if (arb_is_zero(acb_imagref(z)) && arb_is_negative(acb_realref(z)))
+        return 1;
+
+    arf_init(t);
+    mag_init(x);
+    mag_init(y);
+    mag_init(s);
+
+    arf_set_mag(t, arb_radref(acb_realref(z)));
+    arf_add(t, arb_midref(acb_realref(z)), t, MAG_BITS, ARF_RND_CEIL);
+
+    if (arf_sgn(t) >= 0)
     {
-            acb_init(z1);
+        res = 0;
-
-            arb_sqrt_ui(acb_imagref(z1), 3, wp);
-            arb_set_si(acb_realref(z1), -1);
-            acb_mul_2exp_si(z1, z1, -1);
-
-            /* multiply by exp(-2 pi i / 3) in upper half plane
-               and by exp(2 pi i / 3) in lower half plane, to stay close
-               to positive reals */
-            if (arf_sgn(zim) >= 0)
-                acb_conj(z1, z1);
-
-            acb_mul(z1, z1, z, wp);
-            acb_neg(zeta, zeta);   /* same effect regardless of exp(+/-2 pi i/3) */
-
-            acb_hypgeom_airy_bound_9_7_17(B, z1, zeta);
-            mag_mul_2exp_si(B, B, 1);
-            mag_add(B, B, A);
-
-            acb_clear(z1);
     }
+    else
+    {
+        arf_get_mag_lower(x, t);
+        arb_get_mag(y, acb_imagref(z));
+        mag_set_ui_lower(s, 3);
+        mag_sqrt_lower(s, s);
+        mag_mul_lower(s, s, x);
+        res = mag_cmp(s, y) >= 0;
     }
 
-    acb_get_mag(zhi, z);
+    arf_clear(t);
+    mag_clear(x);
+    mag_clear(y);
+    mag_clear(s);
+
+    return res;
+}
+
+void
+acb_hypgeom_airy_bound(mag_t ai, mag_t aip, mag_t bi, mag_t bip, const acb_t z)
+{
+    acb_t zeta;
+    mag_t A, B, D, zlo, zhi;
+    slong wp;
+    int near_zero;
+
+    /* Fast and slightly tighter bounds for real z <= 0 */
+    /* Todo: could implement bounds for z >= 0 too */
+    if (acb_is_real(z) && arb_is_nonpositive(acb_realref(z)))
+    {
+        mag_init(zlo);
+        mag_init(zhi);
+        mag_init(A);
+        mag_init(B);
+        mag_init(D);
+
+        if (ai != NULL || bi != NULL)
+        {
             acb_get_mag_lower(zlo, z);
+            mag_rsqrt(A, zlo);
+            mag_sqrt(A, A);
+            mag_mul_ui(A, A, 150);
+            mag_set_ui(D, 160);
+            mag_min(A, A, D);
+            mag_mul_2exp_si(A, A, -8);
+
+            if (ai != NULL) mag_set(ai, A);
+            if (bi != NULL) mag_set(bi, A);
+        }
+
+        if (aip != NULL || bip != NULL)
+        {
+            acb_get_mag(zhi, z);
+            mag_sqrt(A, zhi);
+            mag_sqrt(A, A);
+            mag_mul_ui(A, A, 150);
+            mag_set_ui(D, 160);
+            mag_max(B, A, D);
+            mag_mul_2exp_si(B, B, -8);
+            mag_set_ui(D, 67);
+            mag_max(A, A, D);
+            mag_mul_2exp_si(A, A, -8);
+
+            if (aip != NULL) mag_set(aip, A);
+            if (bip != NULL) mag_set(bip, B);
+        }
+
+        mag_clear(zlo);
+        mag_clear(zhi);
+        mag_clear(A);
+        mag_clear(B);
+        mag_clear(D);
+        return;
+    }
+
+    acb_init(zeta);
+    mag_init(A);
+    mag_init(B);
+    mag_init(D);
+    mag_init(zlo);
+    mag_init(zhi);
+
+    wp = MAG_BITS * 2;
+
+    acb_get_mag_lower(zlo, z);
+    acb_get_mag(zhi, z);
+
+    if (mag_cmp_2exp_si(zhi, 0) <= 0)
+    {
+        /* inside unit circle -- don't look at asymptotics */
+        if (ai != NULL) mag_set_ui_2exp_si(ai, 159, -8);
+        if (aip != NULL) mag_set_ui_2exp_si(aip, 125, -8);
+        if (bi != NULL) mag_set_ui_2exp_si(bi, 310, -8);
+        if (bip != NULL) mag_set_ui_2exp_si(bip, 239, -8);
+    }
+    else
+    {
+        /* look at asymptotics outside unit circle */
+        near_zero = mag_cmp_2exp_si(zlo, 0) <= 0;
+        if (near_zero)
+            mag_one(zlo);
+
+        if (arg_lt_2pi3_fast(z))
+        {
+            acb_hypgeom_airy_bound_arg_le_2pi3(A, (bi != NULL || bip != NULL) ? B : NULL, z, wp);
+        }
+        else if (arg_gt_2pi3_fast(z))
+        {
+            acb_hypgeom_airy_bound_arg_ge_2pi3(A, (bi != NULL || bip != NULL) ? B : NULL, z, wp);
+        }
+        else
+        {
+            mag_t A2, B2;
+
+            mag_init(A2);
+            mag_init(B2);
+
+            acb_hypgeom_airy_bound_arg_le_2pi3(A, (bi != NULL || bip != NULL) ? B : NULL, z, wp);
+            acb_hypgeom_airy_bound_arg_ge_2pi3(A2, (bi != NULL || bip != NULL) ? B2 : NULL, z, wp);
+
+            mag_max(A, A, A2);
+            mag_max(B, B, A2);
+
+            mag_clear(A2);
+            mag_clear(B2);
+        }
 
         /* bound |z|^(1/4) */
         mag_sqrt(zhi, zhi);
@@ -228,6 +369,35 @@ acb_hypgeom_airy_bound(mag_t ai, mag_t aip, mag_t bi, mag_t bip, const acb_t z)
         if (bi != NULL) mag_mul(bi, B, zlo);
         if (bip != NULL) mag_mul(bip, B, zhi);
 
+        if (near_zero)
+        {
+            /* max by bounds on the unit circle */
+            if (ai != NULL)
+            {
+                mag_set_ui_2exp_si(D, 159, -8);
+                mag_max(ai, ai, D);
+            }
+
+            if (aip != NULL)
+            {
+                mag_set_ui_2exp_si(D, 125, -8);
+                mag_max(aip, aip, D);
+            }
+
+            if (bi != NULL)
+            {
+                mag_set_ui_2exp_si(D, 310, -8);
+                mag_max(bi, bi, D);
+            }
+
+            if (bip != NULL)
+            {
+                mag_set_ui_2exp_si(D, 239, -8);
+                mag_max(bip, bip, D);
+            }
+        }
+    }
+
     acb_clear(zeta);
     mag_clear(A);
     mag_clear(B);

examples/integrals.c

@@ -542,11 +542,25 @@ f_erf_bent(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
     return 0;
 }
 
+/* f(z) = Ai(z) */
+int
+f_airy_ai(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
+{
+    if (order > 1)
+        flint_abort();  /* Would be needed for Taylor method. */
+
+
+    acb_hypgeom_airy(res, NULL, NULL, NULL, z, prec);
+
+    return 0;
+}
+
+
 /* ------------------------------------------------------------------------- */
 /*  Main test program                                                        */
 /* ------------------------------------------------------------------------- */
 
-#define NUM_INTEGRALS 26
+#define NUM_INTEGRALS 27
 
 const char * descr[NUM_INTEGRALS] =
 {
@@ -576,6 +590,7 @@ const char * descr[NUM_INTEGRALS] =
     "int_0^{1000} W_0(x) dx",
     "int_0^pi max(sin(x), cos(x)) dx",
     "int_{-1}^1 erf(x/sqrt(0.0002)*0.5+1.5)*exp(-x) dx",
+    "int_{-10}^10 Ai(x) dx",
 };
 
 int main(int argc, char *argv[])
@@ -975,6 +990,12 @@ int main(int argc, char *argv[])
                 acb_calc_integrate(s, f_erf_bent, NULL, a, b, goal, tol, options, prec);
                 break;
 
+            case 26:
+                acb_set_si(a, -10);
+                acb_set_si(b, 10);
+                acb_calc_integrate(s, f_airy_ai, NULL, a, b, goal, tol, options, prec);
+                break;
+
             default:
                 abort();
             }