From d64055fb3f65f474c0e3d3e8ad1dd58588db9927 Mon Sep 17 00:00:00 2001 From: fredrik Date: Sun, 11 Feb 2018 22:28:49 +0100 Subject: [PATCH] attempt to improve bounds for Airy functions --- acb_hypgeom/airy.c | 42 +++- acb_hypgeom/airy_bound.c | 440 +++++++++++++++++++++++++++------------ examples/integrals.c | 23 +- 3 files changed, 362 insertions(+), 143 deletions(-) diff --git a/acb_hypgeom/airy.c b/acb_hypgeom/airy.c index c0af1148..47246266 100644 --- a/acb_hypgeom/airy.c +++ b/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, - const acb_t z, slong n, slong prec) +acb_hypgeom_airy_prop(acb_t ai, acb_t aip, acb_t bi, acb_t bip, + 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); - acb_hypgeom_airy_direct(ai, aip, bi, bip, zz, n, prec); + 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; } - acb_hypgeom_airy_direct(ai, aip, bi, bip, z, n, wp); + 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); - acb_hypgeom_airy_direct(ai, aip, bi, bip, z, n, wp); + + 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 { diff --git a/acb_hypgeom/airy_bound.c b/acb_hypgeom/airy_bound.c index 48d7dde0..c1853d2d 100644 --- a/acb_hypgeom/airy_bound.c +++ b/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? */ -static int -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; -} - +/* Implements DLMF 9.7.17. We assume |zeta| >= 1/2 and |arg(z)| <= 2 pi/3 here, + ignoring the smaller points which must be dealt with separately. */ 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 */ - 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); - } + /* simplified bound assuming |zeta| >= 1/2 */ + mag_one(D); /* exp(-zeta) / (2 sqrt(pi)) * (1 + D / |zeta|) */ @@ -105,21 +69,238 @@ acb_hypgeom_airy_bound_9_7_17(mag_t bound, const acb_t z, const acb_t zeta) mag_clear(zeta_lower); } +void +acb_hypgeom_airy_bound_arg_le_2pi3(mag_t A, mag_t B, const acb_t z, slong wp) +{ + acb_t zeta, z1; + + acb_init(zeta); + acb_init(z1); + + acb_set_round(zeta, z, wp); + acb_sqrt(zeta, zeta, wp); + acb_cube(zeta, zeta, wp); + acb_mul_2exp_si(zeta, zeta, 1); + acb_div_ui(zeta, zeta, 3, wp); + + acb_hypgeom_airy_bound_9_7_17(A, z, zeta); + + /* 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); + + arb_sqrt_ui(acb_imagref(z1), 3, wp); + arb_one(acb_realref(z1)); + acb_mul_2exp_si(z1, z1, -1); + acb_conj(z2, z1); + + acb_neg_round(zeta, z, wp); + acb_mul(z1, z1, zeta, wp); + + acb_sqrt(zeta, zeta, wp); + acb_cube(zeta, zeta, wp); + acb_mul_2exp_si(zeta, zeta, 1); + acb_div_ui(zeta, zeta, 3, wp); + acb_mul_onei(zeta, zeta); + + acb_hypgeom_airy_bound_9_7_17(A, z1, zeta); + + /* conjugate symmetry */ + if (arb_is_zero(acb_imagref(z))) + { + mag_mul_2exp_si(A, A, 1); + } + 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 + { + arf_get_mag(x, t); + arb_get_mag_lower(y, acb_imagref(z)); + mag_set_ui(s, 3); + mag_sqrt(s, s); + mag_mul(s, s, x); + res = mag_cmp(s, y) <= 0; + } + + arf_clear(t); + mag_clear(x); + mag_clear(y); + mag_clear(s); + + return res; +} + +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) + { + res = 0; + } + 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; + } + + 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; - acb_t z1, z2; - arf_srcptr zre, zim; mag_t A, B, D, zlo, zhi; slong wp; + int near_zero; - if (acb_contains_zero(z)) + /* 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))) { - if (ai != NULL) mag_inf(ai); - if (aip != NULL) mag_inf(aip); - if (bi != NULL) mag_inf(bi); - if (bip != NULL) mag_inf(bip); + 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; } @@ -131,103 +312,92 @@ acb_hypgeom_airy_bound(mag_t ai, mag_t aip, mag_t bi, mag_t bip, const acb_t z) mag_init(zhi); wp = MAG_BITS * 2; - zre = arb_midref(acb_realref(z)); - zim = arb_midref(acb_imagref(z)); - /* near the negative half line, use - |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) + acb_get_mag_lower(zlo, z); + acb_get_mag(zhi, z); + + if (mag_cmp_2exp_si(zhi, 0) <= 0) { - acb_init(z1); - acb_init(z2); - - arb_sqrt_ui(acb_imagref(z1), 3, wp); - arb_one(acb_realref(z1)); - acb_mul_2exp_si(z1, z1, -1); - acb_conj(z2, z1); - - acb_neg_round(zeta, z, wp); - acb_mul(z1, z1, zeta, wp); - - acb_sqrt(zeta, zeta, wp); - acb_cube(zeta, zeta, wp); - acb_mul_2exp_si(zeta, zeta, 1); - acb_div_ui(zeta, zeta, 3, wp); - acb_mul_onei(zeta, zeta); - - acb_hypgeom_airy_bound_9_7_17(A, z1, zeta); - - /* conjugate symmetry */ - if (arb_is_zero(acb_imagref(z))) - { - mag_mul_2exp_si(A, A, 1); - } - else - { - 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_set(B, A); - - acb_clear(z1); - acb_clear(z2); + /* 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 { - acb_set_round(zeta, z, wp); - acb_sqrt(zeta, zeta, wp); - acb_cube(zeta, zeta, wp); - acb_mul_2exp_si(zeta, zeta, 1); - acb_div_ui(zeta, zeta, 3, wp); + /* look at asymptotics outside unit circle */ + near_zero = mag_cmp_2exp_si(zlo, 0) <= 0; + if (near_zero) + mag_one(zlo); - acb_hypgeom_airy_bound_9_7_17(A, z, zeta); - - /* 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 (bi != NULL || bip != NULL) + if (arg_lt_2pi3_fast(z)) { - acb_init(z1); + 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; - arb_sqrt_ui(acb_imagref(z1), 3, wp); - arb_set_si(acb_realref(z1), -1); - acb_mul_2exp_si(z1, z1, -1); + mag_init(A2); + mag_init(B2); - /* 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_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); - acb_mul(z1, z1, z, wp); - acb_neg(zeta, zeta); /* same effect regardless of exp(+/-2 pi i/3) */ + mag_max(A, A, A2); + mag_max(B, B, A2); - acb_hypgeom_airy_bound_9_7_17(B, z1, zeta); - mag_mul_2exp_si(B, B, 1); - mag_add(B, B, A); + mag_clear(A2); + mag_clear(B2); + } - acb_clear(z1); + /* bound |z|^(1/4) */ + mag_sqrt(zhi, zhi); + mag_sqrt(zhi, zhi); + + /* bound |z|^(-1/4) */ + mag_rsqrt(zlo, zlo); + mag_sqrt(zlo, zlo); + + if (ai != NULL) mag_mul(ai, A, zlo); + if (aip != NULL) mag_mul(aip, A, zhi); + 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_get_mag(zhi, z); - acb_get_mag_lower(zlo, z); - - /* bound |z|^(1/4) */ - mag_sqrt(zhi, zhi); - mag_sqrt(zhi, zhi); - - /* bound |z|^(-1/4) */ - mag_rsqrt(zlo, zlo); - mag_sqrt(zlo, zlo); - - if (ai != NULL) mag_mul(ai, A, zlo); - if (aip != NULL) mag_mul(aip, A, zhi); - if (bi != NULL) mag_mul(bi, B, zlo); - if (bip != NULL) mag_mul(bip, B, zhi); - acb_clear(zeta); mag_clear(A); mag_clear(B); diff --git a/examples/integrals.c b/examples/integrals.c index 63e8c448..e303d3d0 100644 --- a/examples/integrals.c +++ 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(); }