From 9eed54684cf56461d787f21731a0f943d5f528a3 Mon Sep 17 00:00:00 2001 From: Fredrik Johansson Date: Thu, 19 Nov 2015 11:10:22 +0100 Subject: [PATCH] initial code for Airy functions --- acb_hypgeom.h | 5 + acb_hypgeom/airy.c | 294 ++++++++++++++++ acb_hypgeom/airy_asymp.c | 601 ++++++++++++++++++++++++++++++++ acb_hypgeom/airy_bound.c | 252 +++++++++++++ acb_hypgeom/airy_direct.c | 301 ++++++++++++++++ acb_hypgeom/test/t-airy.c | 173 +++++++++ acb_hypgeom/test/t-airy_bound.c | 158 +++++++++ doc/source/acb_hypgeom.rst | 142 ++++++++ 8 files changed, 1926 insertions(+) create mode 100644 acb_hypgeom/airy.c create mode 100644 acb_hypgeom/airy_asymp.c create mode 100644 acb_hypgeom/airy_bound.c create mode 100644 acb_hypgeom/airy_direct.c create mode 100644 acb_hypgeom/test/t-airy.c create mode 100644 acb_hypgeom/test/t-airy_bound.c diff --git a/acb_hypgeom.h b/acb_hypgeom.h index 1554f6d8..f9ab2de4 100644 --- a/acb_hypgeom.h +++ b/acb_hypgeom.h @@ -113,6 +113,11 @@ void acb_hypgeom_0f1_asymp(acb_t res, const acb_t a, const acb_t z, int regulari void acb_hypgeom_0f1_direct(acb_t res, const acb_t a, const acb_t z, int regularized, slong prec); void acb_hypgeom_0f1(acb_t res, const acb_t a, const acb_t z, int regularized, slong prec); +void acb_hypgeom_airy_bound(mag_t ai, mag_t aip, mag_t bi, mag_t bip, const acb_t z); +void acb_hypgeom_airy_asymp(acb_t ai, acb_t aip, acb_t bi, acb_t bip, const acb_t z, slong n, slong prec); +void acb_hypgeom_airy_direct(acb_t ai, acb_t aip, acb_t bi, acb_t bip, const acb_t z, slong n, slong prec); +void acb_hypgeom_airy(acb_t ai, acb_t aip, acb_t bi, acb_t bip, const acb_t z, slong prec); + void acb_hypgeom_gamma_upper_asymp(acb_t res, const acb_t s, const acb_t z, int modified, slong prec); void acb_hypgeom_gamma_upper_1f1a(acb_t res, const acb_t s, const acb_t z, int modified, slong prec); void acb_hypgeom_gamma_upper_1f1b(acb_t res, const acb_t s, const acb_t z, int modified, slong prec); diff --git a/acb_hypgeom/airy.c b/acb_hypgeom/airy.c new file mode 100644 index 00000000..ea601fc4 --- /dev/null +++ b/acb_hypgeom/airy.c @@ -0,0 +1,294 @@ +/*============================================================================= + + This file is part of ARB. + + ARB is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + ARB is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ARB; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2015 Fredrik Johansson + +******************************************************************************/ + +#include "acb_hypgeom.h" +#include "double_extras.h" + +#define LOG2 0.69314718055994530942 +#define EXP1 2.7182818284590452354 + +static const double small_log_tab[] = { + 0.0, 0.0, 0.693147180559945309, + 1.09861228866810969, 1.38629436111989062, 1.60943791243410037, + 1.791759469228055, 1.94591014905531331, 2.07944154167983593, + 2.19722457733621938, 2.30258509299404568, 2.39789527279837054, + 2.48490664978800031, 2.56494935746153674, 2.63905732961525861, + 2.70805020110221007, 2.77258872223978124, 2.83321334405621608, + 2.89037175789616469, 2.94443897916644046, 2.99573227355399099, + 3.044522437723423, 3.09104245335831585, 3.13549421592914969, + 3.17805383034794562, 3.21887582486820075, 3.25809653802148205, + 3.29583686600432907, 3.33220451017520392, 3.36729582998647403, + 3.40119738166215538, 3.43398720448514625, 3.46573590279972655, + 3.49650756146648024, 3.52636052461616139, 3.55534806148941368, + 3.58351893845611, 3.61091791264422444, 3.63758615972638577, + 3.66356164612964643, 3.6888794541139363, 3.7135720667043078, + 3.73766961828336831, 3.76120011569356242, 3.78418963391826116, + 3.80666248977031976, 3.828641396489095, 3.85014760171005859, + 3.87120101090789093, 3.89182029811062661, 3.91202300542814606, + 3.93182563272432577, 3.95124371858142735, 3.97029191355212183, + 3.98898404656427438, 4.00733318523247092, 4.02535169073514923, + 4.04305126783455015, 4.06044301054641934, 4.07753744390571945, + 4.09434456222210068, 4.11087386417331125, 4.12713438504509156, + 4.14313472639153269, +}; + +static slong +asymp_pick_terms(double prec, double logz) +{ + double logk, log2term, log2term_prev; + slong k; + + log2term_prev = 0.0; + + for (k = 1; ; k++) + { + logk = k < 64 ? small_log_tab[k] : log(k); + + log2term = -1.3257480647361593990 - 0.72134752044448170368*logk + + k * (-1.8577325401678072259 + 1.4426950408889634074*logk + - 2.1640425613334451110*logz); + + if (log2term > log2term_prev - 0.5) + return -1; + + if (log2term < -prec) + return k; + } +} + +/* +Accurate estimate of log2(|Ai(z)|) or log2(|Bi(z)|), given z = x + yi. +Should subtract 0.25*log2(|z| + log2(2*sqrt(pi)) from the output. +*/ +static double +estimate_airy(double x, double y, int ai) +{ + double r, t, sgn; + + r = x; + t = x * ((x * x) - 3.0 * (y * y)); + y = y * (3.0 * (x * x) - (y * y)); + x = t * (1.0 / 9.0); + y = y * (1.0 / 9.0); + + sgn = 1.0; + + if (r > 0.0 && x > 0.0) + { + t = sqrt(x * x + y * y) + x; + + if (ai) sgn = -1.0; + } + else + { + x = -x; + + if (x > 0.0 && x > 1e6 * fabs(y)) + t = y * y / (2.0 * x); + else + t = sqrt(x * x + y * y) - x; + } + + return sgn * sqrt(0.5 * t) * 2.8853900817779268147; +} + +/* 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, long prec) +{ + mag_t aib, aipb, bib, bipb, zb, rad; + acb_t zz; + int real; + + mag_init(aib); + mag_init(aipb); + mag_init(bib); + mag_init(bipb); + mag_init(zb); + mag_init(rad); + acb_init(zz); + + real = acb_is_real(z); + arf_set(arb_midref(acb_realref(zz)), arb_midref(acb_realref(z))); + arf_set(arb_midref(acb_imagref(zz)), arb_midref(acb_imagref(z))); + mag_hypot(rad, arb_radref(acb_realref(z)), arb_radref(acb_imagref(z))); + 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 (ai != NULL) + { + mag_mul(aipb, aipb, rad); + if (real) + arb_add_error_mag(acb_realref(ai), aipb); + else + acb_add_error_mag(ai, aipb); + } + + if (aip != NULL) + { + mag_mul(aib, aib, rad); + mag_mul(aib, aib, zb); /* |Ai''(z)| = |z Ai(z)| */ + if (real) + arb_add_error_mag(acb_realref(aip), aib); + else + acb_add_error_mag(aip, aib); + } + + if (bi != NULL) + { + mag_mul(bipb, bipb, rad); + if (real) + arb_add_error_mag(acb_realref(bi), bipb); + else + acb_add_error_mag(bi, bipb); + } + + if (bip != NULL) + { + mag_mul(bib, bib, rad); + mag_mul(bib, bib, zb); /* |Bi''(z)| = |z Bi(z)| */ + if (real) + arb_add_error_mag(acb_realref(bip), bib); + else + acb_add_error_mag(bip, bib); + } + + mag_clear(aib); + mag_clear(aipb); + mag_clear(bib); + mag_clear(bipb); + mag_clear(zb); + mag_clear(rad); + acb_clear(zz); +} + +void +acb_hypgeom_airy(acb_t ai, acb_t aip, acb_t bi, acb_t bip, const acb_t z, long prec) +{ + arf_srcptr re, im; + double x, y, t, zmag, z15, term_est, airy_est, abstol; + slong n, wp; + + if (!acb_is_finite(z)) + { + if (ai != NULL) acb_indeterminate(ai); + if (aip != NULL) acb_indeterminate(aip); + if (bi != NULL) acb_indeterminate(bi); + if (bip != NULL) acb_indeterminate(bip); + return; + } + + re = arb_midref(acb_realref(z)); + im = arb_midref(acb_imagref(z)); + wp = prec * 1.03 + 15; + + /* tiny input -- use direct method and pick n without underflowing */ + if (arf_cmpabs_2exp_si(re, -64) < 0 && arf_cmpabs_2exp_si(im, -64) < 0) + { + if (arf_cmpabs_2exp_si(re, -wp) < 0 && arf_cmpabs_2exp_si(im, -wp) < 0) + { + n = 1; /* very tiny input */ + } + else + { + if (arf_cmpabs(re, im) > 0) + zmag = fmpz_get_d(ARF_EXPREF(re)); + else + zmag = fmpz_get_d(ARF_EXPREF(im)); + zmag = (zmag + 1) * (1.0 / LOG2); + n = wp / (-zmag) + 1; + } + + acb_hypgeom_airy_direct(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)) + { + if (arf_cmpabs_2exp_si(re, prec) > 0 || arf_cmpabs_2exp_si(im, prec) > 0) + { + n = 1; /* very huge input */ + } + else + { + x = fmpz_get_d(ARF_EXPREF(re)); + y = fmpz_get_d(ARF_EXPREF(im)); + zmag = (FLINT_MAX(x, y) - 2) * (1.0 / LOG2); + n = asymp_pick_terms(wp, zmag); + n = FLINT_MAX(n, 1); + } + + acb_hypgeom_airy_asymp(ai, aip, bi, bip, z, n, wp); + } + else /* moderate input */ + { + x = arf_get_d(re, ARF_RND_DOWN); + y = arf_get_d(im, ARF_RND_DOWN); + + zmag = sqrt(x * x + y * y); + z15 = zmag * sqrt(zmag); + + if (zmag >= 4.0 && (n = asymp_pick_terms(wp, log(zmag))) != -1) + { + acb_hypgeom_airy_asymp(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); + } + else + { + /* estimate largest term: log2(exp(2(z^3/9)^(1/2))) */ + term_est = 0.96179669392597560491 * z15; + + /* estimate the smaller of Ai and Bi */ + airy_est = estimate_airy(x, y, (ai != NULL || aip != NULL)); + + /* estimate absolute tolerance and necessary working precision */ + abstol = airy_est - wp; + wp = wp + term_est - airy_est; + wp = FLINT_MAX(wp, 10); + + t = 3 * (-abstol * LOG2) / (2 * z15 * EXP1); + t = (-abstol * 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); + } + } + + if (ai != NULL) acb_set_round(ai, ai, prec); + if (aip != NULL) acb_set_round(aip, aip, prec); + if (bi != NULL) acb_set_round(bi, bi, prec); + if (bip != NULL) acb_set_round(bip, bip, prec); +} + diff --git a/acb_hypgeom/airy_asymp.c b/acb_hypgeom/airy_asymp.c new file mode 100644 index 00000000..37769e62 --- /dev/null +++ b/acb_hypgeom/airy_asymp.c @@ -0,0 +1,601 @@ +/*============================================================================= + + This file is part of ARB. + + ARB is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + ARB is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ARB; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2015 Fredrik Johansson + +******************************************************************************/ + +#include "acb_hypgeom.h" + +void acb_hypgeom_mag_chi(mag_t chi, ulong n); + +static int +arg_lt_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_positive(acb_imagref(zeta))) + return 1; + + if (arb_is_negative(acb_imagref(z)) && + arb_is_negative(acb_imagref(zeta))) + return 1; + + return 0; +} + +/* assuming a >= b >= c >= d, e >= f */ +void +_acb_mul4div2_ui(acb_t x, ulong a, ulong b, ulong c, ulong d, ulong e, ulong f, slong prec) +{ + if (a < (UWORD(1) << (FLINT_BITS / 4))) + { + acb_mul_ui(x, x, a * b * c * d, prec); + } + else if (a < (UWORD(1) << (FLINT_BITS / 2))) + { + acb_mul_ui(x, x, a * b, prec); + acb_mul_ui(x, x, c * d, prec); + } + else + { + acb_mul_ui(x, x, a, prec); + acb_mul_ui(x, x, b, prec); + acb_mul_ui(x, x, c, prec); + acb_mul_ui(x, x, d, prec); + } + + if (e < (UWORD(1) << (FLINT_BITS / 2))) + { + acb_div_ui(x, x, e * f, prec); + } + else + { + acb_div_ui(x, x, e, prec); + acb_div_ui(x, x, f, prec); + } +} + +void +acb_hypgeom_airy_asymp_sum(acb_t s0even, acb_t s0odd, + acb_t s1even, acb_t s1odd, + mag_t t0n, mag_t t1n, + const acb_t z, const acb_t z2, slong n, slong prec) +{ + slong m, k, j; + acb_ptr z2pow; + + acb_get_mag(t0n, z); + mag_mul_ui(t0n, t0n, 72); + mag_pow_ui(t0n, t0n, n); + mag_one(t1n); + + for (k = 1; k <= n; k++) + { + mag_mul_ui(t0n, t0n, 6 * k - 1); + mag_mul_ui(t0n, t0n, 6 * k - 5); + mag_mul_ui_lower(t1n, t1n, 72 * k); + } + + mag_div(t0n, t0n, t1n); + mag_mul_ui(t1n, t0n, 6 * n + 1); + mag_div_ui(t1n, t1n, 6 * n - 1); + + m = n_sqrt(n / 2); + m = FLINT_MAX(m, 1); + + z2pow = _acb_vec_init(m + 1); + _acb_vec_set_powers(z2pow, z2, m + 1, prec); + + if (s0even != NULL) + { + acb_zero(s0even); + + for (k = n + (n % 2); k >= 0; k -= 2) + { + j = (k / 2) % m; + + if (k < n) + acb_add(s0even, s0even, z2pow + j, prec); + + if (k > 0) + { + _acb_mul4div2_ui(s0even, 6*k-1, 6*k-5, 6*k-7, 6*k-11, k, k-1, prec); + if (j == 0 && k < n) + acb_mul(s0even, s0even, z2pow + m, prec); + } + } + } + + if (s0odd != NULL) + { + acb_zero(s0odd); + + for (k = n + 1 + (n % 2); k >= 1; k -= 2) + { + j = ((k - 1) / 2) % m; + + if (k < n) + acb_add(s0odd, s0odd, z2pow + j, prec); + + if (k > 1) + { + _acb_mul4div2_ui(s0odd, 6*k-1, 6*k-5, 6*k-7, 6*k-11, k, k-1, prec); + if (j == 0 && k < n) + acb_mul(s0odd, s0odd, z2pow + m, prec); + } + else + { + acb_mul(s0odd, s0odd, z, prec); + acb_mul_ui(s0odd, s0odd, 5, prec); + } + } + } + + if (s1even != NULL) + { + acb_zero(s1even); + + for (k = n + (n % 2); k >= 0; k -= 2) + { + j = (k / 2) % m; + + if (k < n) + acb_add(s1even, s1even, z2pow + j, prec); + + if (k > 0) + { + _acb_mul4div2_ui(s1even, 6*k+1, 6*k-5, 6*k-7, FLINT_ABS(6*k-13), k, k-1, prec); + if (k == 2) + acb_neg(s1even, s1even); + if (j == 0 && k < n) + acb_mul(s1even, s1even, z2pow + m, prec); + } + } + } + + if (s1odd != NULL) + { + acb_zero(s1odd); + + for (k = n + 1 + (n % 2); k >= 1; k -= 2) + { + j = ((k - 1) / 2) % m; + + if (k < n) + acb_add(s1odd, s1odd, z2pow + j, prec); + + if (k > 1) + { + _acb_mul4div2_ui(s1odd, 6*k+1, 6*k-5, 6*k-7, 6*k-13, k, k-1, prec); + if (j == 0 && k < n) + acb_mul(s1odd, s1odd, z2pow + m, prec); + } + else + { + acb_mul(s1odd, s1odd, z, prec); + acb_mul_si(s1odd, s1odd, -7, prec); + } + } + } + + _acb_vec_clear(z2pow, m + 1); +} + +void +acb_hypgeom_airy_asymp_bound_factor(mag_t bound, const acb_t z, const acb_t zeta, ulong n) +{ + mag_t t, u, v; + + mag_init(t); + mag_init(u); + mag_init(v); + + if (arb_is_nonnegative(acb_realref(z)) && arb_is_nonnegative(acb_realref(zeta))) + { + /* 2 exp(7 / (36 |zeta|)) */ + mag_set_ui_2exp_si(t, 50, -8); /* bound for 7/36 */ + acb_get_mag_lower(u, zeta); + mag_div(t, t, u); + mag_exp(t, t); + mag_mul_2exp_si(bound, t, 1); + } + else + { + /* 2 exp(7 pi / (72 |zeta|)) */ + mag_set_ui_2exp_si(t, 79, -8); /* bound for 7 pi/72 */ + acb_get_mag_lower(u, zeta); + mag_div(t, t, u); + mag_exp(t, t); + mag_mul_2exp_si(bound, t, 1); + + /* 4 exp(7 pi / (36 |re(zeta)|)) / |cos(arg(zeta))|^n */ + if (!arg_lt_2pi3(z, zeta)) + { + arb_get_mag_lower(u, acb_realref(zeta)); + arb_get_mag(v, acb_imagref(zeta)); + + /* 4 exp(7 pi / (36 u)) < exp(157/(256 u)) */ + mag_set_ui_2exp_si(t, 157, -8); + mag_div(t, t, u); + mag_exp(t, t); + mag_mul_2exp_si(t, t, 2); + + /* |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_pow_ui(v, v, n); + + mag_mul(t, t, v); + mag_max(bound, bound, t); + } + + /* chi(n) */ + acb_hypgeom_mag_chi(t, n); + mag_mul(bound, bound, t); + } + + mag_clear(t); + mag_clear(u); + mag_clear(v); +} + +void +acb_hypgeom_airy_asymp(acb_t ai, acb_t aip, acb_t bi, acb_t bip, const acb_t z, slong n, slong prec) +{ + acb_t t, u, w, z_root, zeta; + acb_t s0even, s0odd, s1even, s1odd, E1, E2; + mag_t err1, err2, erru, errv, errtmp; + int want_d0, want_d1, is_real, upper; + + acb_init(t); + acb_init(u); + acb_init(w); + acb_init(z_root); + acb_init(zeta); + acb_init(s0even); + acb_init(s0odd); + acb_init(s1even); + acb_init(s1odd); + acb_init(E1); + acb_init(E2); + + mag_init(err1); + mag_init(err2); + mag_init(errtmp); + mag_init(erru); + mag_init(errv); + + is_real = acb_is_real(z); + want_d0 = (ai != NULL) || (bi != NULL); + want_d1 = (aip != NULL) || (bip != NULL); + + /* required for some of the error bounds to be valid */ + n = FLINT_MAX(n, 1); + + /* this is just a heuristic; the sectors are validated later */ + if (arf_sgn(arb_midref(acb_realref(z))) >= 0) + { + /* z_root = z^(1/4), zeta = 2 z^(3/2) / 3 */ + acb_sqrt(z_root, z, prec); + acb_cube(zeta, z_root, prec); + acb_sqrt(z_root, z_root, prec); + acb_mul_2exp_si(zeta, zeta, 1); + acb_div_ui(zeta, zeta, 3, prec); + + /* compute bound factor for z = t = z1, zeta = u = zeta1 */ + acb_hypgeom_airy_asymp_bound_factor(err1, z, zeta, n); + + /* this is just a heuristic; the sectors are validated later */ + upper = arf_sgn(arb_midref(acb_imagref(z))) >= 0; + + if (upper) + { + /* check -pi/3 < arg(z) < pi */ + if (arb_is_positive(acb_imagref(z)) || + (arb_is_positive(acb_realref(z)) && + arb_is_positive(acb_realref(zeta)))) + { + arb_one(acb_realref(w)); + arb_sqrt_ui(acb_imagref(w), 3, prec); + acb_mul_2exp_si(w, w, -1); + acb_neg(w, w); + acb_mul(t, w, z, prec); + acb_neg(u, zeta); + acb_hypgeom_airy_asymp_bound_factor(err2, t, u, n); + } + else + { + mag_inf(err2); + } + } + else + { + /* check -pi < arg(z) < pi/3 */ + if (arb_is_negative(acb_imagref(z)) || + (arb_is_positive(acb_realref(z)) && + arb_is_positive(acb_realref(zeta)))) + { + arb_set_si(acb_realref(w), -1); + arb_sqrt_ui(acb_imagref(w), 3, prec); + acb_mul_2exp_si(w, w, -1); + acb_mul(t, w, z, prec); + acb_neg(u, zeta); + acb_hypgeom_airy_asymp_bound_factor(err2, t, u, n); + } + else + { + mag_inf(err2); + } + } + + acb_mul_ui(t, zeta, 72, prec); + acb_inv(t, t, prec); + acb_mul(u, t, t, prec); + + acb_hypgeom_airy_asymp_sum(want_d0 ? s0even : NULL, + want_d0 ? s0odd : NULL, + want_d1 ? s1even : NULL, + want_d1 ? s1odd : NULL, + erru, errv, t, u, n, prec); + + /* exponentials */ + acb_exp_invexp(E2, E1, zeta, prec); + + /* common prefactor */ + acb_const_pi(w, prec); + acb_sqrt(w, w, prec); + acb_mul_2exp_si(w, w, 1); + + if (ai != NULL || bi != NULL) + { + /* t = sum (-1)^k u(k) / (zeta)^k + error */ + acb_sub(t, s0even, s0odd, prec); + mag_mul(errtmp, err1, erru); + acb_add_error_mag(t, errtmp); + + /* u = sum (-1)^k u(k) / (-zeta)^k = sum u(k) / (zeta)^k + error */ + acb_add(u, s0even, s0odd, prec); + mag_mul(errtmp, err2, erru); + acb_add_error_mag(u, errtmp); + + if (ai != NULL) + { + acb_mul(ai, t, E1, prec); + } + + if (bi != NULL) + { + acb_mul(t, t, E1, prec); + acb_mul(u, u, E2, prec); + acb_mul_2exp_si(u, u, 1); + + if (upper) + acb_mul_onei(t, t); + else + acb_div_onei(t, t); + + acb_add(bi, t, u, prec); + } + + acb_mul(t, w, z_root, prec); + acb_inv(t, t, prec); + + if (ai != NULL) acb_mul(ai, ai, t, prec); + if (bi != NULL) acb_mul(bi, bi, t, prec); + } + + if (aip != NULL || bip != NULL) + { + /* t = sum (-1)^k v(k) / (zeta)^k + error */ + acb_sub(t, s1even, s1odd, prec); + mag_mul(errtmp, err1, errv); + acb_add_error_mag(t, errtmp); + + /* u = sum (-1)^k v(k) / (-zeta)^k = sum v(k) / (zeta)^k + error */ + acb_add(u, s1even, s1odd, prec); + mag_mul(errtmp, err2, errv); + acb_add_error_mag(u, errtmp); + + if (aip != NULL) + { + acb_mul(aip, t, E1, prec); + acb_neg(aip, aip); + } + + if (bip != NULL) + { + acb_mul(t, t, E1, prec); + acb_mul(u, u, E2, prec); + acb_mul_2exp_si(u, u, 1); + + if (upper) + acb_div_onei(t, t); + else + acb_mul_onei(t, t); + + acb_add(bip, t, u, prec); + } + + acb_inv(t, w, prec); + acb_mul(t, t, z_root, prec); + + if (aip != NULL) acb_mul(aip, aip, t, prec); + if (bip != NULL) acb_mul(bip, bip, t, prec); + } + } + else + { + /* z_root = (-z)^(1/4), zeta = 2 (-z)^(3/2) / 3 */ + acb_neg(t, z); + acb_sqrt(z_root, t, prec); + acb_cube(zeta, z_root, prec); + acb_sqrt(z_root, z_root, prec); + acb_mul_2exp_si(zeta, zeta, 1); + acb_div_ui(zeta, zeta, 3, prec); + + if (arg_lt_2pi3(t, zeta)) + { + /* compute bound factor for i zeta */ + arb_one(acb_realref(w)); + arb_sqrt_ui(acb_imagref(w), 3, prec); + acb_mul_2exp_si(w, w, -1); + acb_mul(t, z, w, prec); + acb_neg(t, t); + acb_mul_onei(u, zeta); + acb_hypgeom_airy_asymp_bound_factor(err1, t, u, n); + + /* compute bound factor for -i zeta */ + acb_conj(w, w); + acb_mul(t, z, w, prec); + acb_neg(t, t); + acb_div_onei(u, zeta); + acb_hypgeom_airy_asymp_bound_factor(err2, t, u, n); + } + else + { + mag_inf(err1); + mag_inf(err2); + } + + /* t = 1/(72 i zeta), u = (1/(72 i zeta))^2 */ + acb_mul_onei(t, zeta); + acb_mul_ui(t, t, 72, prec); + acb_inv(t, t, prec); + acb_mul(u, t, t, prec); + + acb_hypgeom_airy_asymp_sum(want_d0 ? s0even : NULL, + want_d0 ? s0odd : NULL, + want_d1 ? s1even : NULL, + want_d1 ? s1odd : NULL, + erru, errv, t, u, n, prec); + + /* exponentials */ + acb_const_pi(t, prec); + acb_mul_2exp_si(t, t, -2); + acb_sub(t, zeta, t, prec); + acb_mul_onei(t, t); + acb_exp_invexp(E2, E1, t, prec); + + /* common prefactor */ + acb_const_pi(w, prec); + acb_sqrt(w, w, prec); + acb_mul_2exp_si(w, w, 1); + + if (ai != NULL || bi != NULL) + { + /* t = sum (-1)^k u(k) / (i zeta)^k + error */ + acb_sub(t, s0even, s0odd, prec); + mag_mul(errtmp, err1, erru); + acb_add_error_mag(t, errtmp); + + /* w = sum (-1)^k u(k) / (-i zeta)^k = sum u(k) / (i zeta)^k + error */ + acb_add(u, s0even, s0odd, prec); + mag_mul(errtmp, err2, erru); + acb_add_error_mag(u, errtmp); + + if (ai != NULL) + { + acb_mul(ai, t, E1, prec); + acb_addmul(ai, u, E2, prec); + } + + if (bi != NULL) + { + /* (-i E1) t + (+i E2) u */ + acb_mul(bi, u, E2, prec); + acb_submul(bi, t, E1, prec); + acb_mul_onei(bi, bi); + } + + acb_mul(t, w, z_root, prec); + acb_inv(t, t, prec); + + if (ai != NULL) acb_mul(ai, ai, t, prec); + if (bi != NULL) acb_mul(bi, bi, t, prec); + } + + if (aip != NULL || bip != NULL) + { + /* t = sum (-1)^k v(k) / (i zeta)^k + error */ + acb_sub(t, s1even, s1odd, prec); + mag_mul(errtmp, err1, errv); + acb_add_error_mag(t, errtmp); + + /* u = sum (-1)^k v(k) / (-i zeta)^k = sum v(k) / (i zeta)^k + error */ + acb_add(u, s1even, s1odd, prec); + mag_mul(errtmp, err2, errv); + acb_add_error_mag(u, errtmp); + + if (bip != NULL) + { + acb_mul(bip, t, E1, prec); + acb_addmul(bip, u, E2, prec); + } + + if (aip != NULL) + { + /* i E1 t - i E2 u */ + acb_mul(aip, t, E1, prec); + acb_submul(aip, u, E2, prec); + acb_mul_onei(aip, aip); + } + + acb_inv(t, w, prec); + acb_mul(t, t, z_root, prec); + + if (aip != NULL) acb_mul(aip, aip, t, prec); + if (bip != NULL) acb_mul(bip, bip, t, prec); + } + } + + if (is_real) + { + if (ai != NULL) arb_zero(acb_imagref(ai)); + if (aip != NULL) arb_zero(acb_imagref(aip)); + if (bi != NULL) arb_zero(acb_imagref(bi)); + if (bip != NULL) arb_zero(acb_imagref(bip)); + } + + acb_clear(t); + acb_clear(u); + acb_clear(w); + acb_clear(z_root); + acb_clear(zeta); + acb_clear(s0even); + acb_clear(s0odd); + acb_clear(s1even); + acb_clear(s1odd); + acb_clear(E1); + acb_clear(E2); + + mag_clear(err1); + mag_clear(err2); + mag_clear(errtmp); + mag_clear(erru); + mag_clear(errv); +} + diff --git a/acb_hypgeom/airy_bound.c b/acb_hypgeom/airy_bound.c new file mode 100644 index 00000000..127a488b --- /dev/null +++ b/acb_hypgeom/airy_bound.c @@ -0,0 +1,252 @@ +/*============================================================================= + + This file is part of ARB. + + ARB is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + ARB is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ARB; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2015 Fredrik Johansson + +******************************************************************************/ + +#include "acb_hypgeom.h" + +void +arb_bound_exp_neg(mag_t b, const arb_t x) +{ + arb_t t; + arb_init(t); + arf_set_mag(arb_midref(t), arb_radref(x)); + arf_sub(arb_midref(t), arb_midref(x), arb_midref(t), MAG_BITS, ARF_RND_FLOOR); + arf_neg(arb_midref(t), arb_midref(t)); + arb_exp(t, t, MAG_BITS); + arb_get_mag(b, t); + 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; +} + +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_init(D); + mag_init(t); + mag_init(u); + mag_init(v); + mag_init(zeta_lower); + + acb_get_mag_lower(zeta_lower, zeta); + + /* 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); + } + + /* exp(-zeta) / (2 sqrt(pi)) * (1 + D / |zeta|) */ + + /* t = exp(-zeta) / (2 sqrt(pi)) < exp(-zeta) * (73/256) */ + arb_bound_exp_neg(t, acb_realref(zeta)); + mag_mul_ui(t, t, 73); + mag_mul_2exp_si(t, t, -8); + + /* u = 1 + D / |zeta| */ + mag_div(u, D, zeta_lower); + mag_one(v); + mag_add(u, u, v); + + mag_mul(bound, t, u); + + mag_clear(D); + mag_clear(t); + mag_clear(u); + mag_clear(v); + mag_clear(zeta_lower); +} + +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; + + 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); + mag_init(B); + mag_init(D); + mag_init(zlo); + 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_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); + } + 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); + + 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) + { + acb_init(z1); + + 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); + } + } + + 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); + mag_clear(D); + mag_clear(zlo); + mag_clear(zhi); +} + diff --git a/acb_hypgeom/airy_direct.c b/acb_hypgeom/airy_direct.c new file mode 100644 index 00000000..855fc318 --- /dev/null +++ b/acb_hypgeom/airy_direct.c @@ -0,0 +1,301 @@ +/*============================================================================= + + This file is part of ARB. + + ARB is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + ARB is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ARB; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2015 Fredrik Johansson + +******************************************************************************/ + +#include "acb_hypgeom.h" + +void +arb_const_airy_ai0_eval(arb_t y, slong prec) +{ + arb_t t; fmpq_t v; arb_init(t); fmpq_init(v); + arb_set_ui(y, 3); + arb_root(y, y, 3, prec + 5); arb_mul(y, y, y, prec + 5); + fmpq_set_si(v, 2, 3); arb_gamma_fmpq(t, v, prec + 5); + arb_mul(y, y, t, prec + 5); + arb_inv(y, y, prec); + arb_clear(t); fmpq_clear(v); +} + +void +arb_const_airy_ai1_eval(arb_t y, slong prec) +{ + arb_t t; fmpq_t v; arb_init(t); fmpq_init(v); + arb_set_ui(y, 3); + arb_root(y, y, 3, prec + 5); + fmpq_set_si(v, 1, 3); arb_gamma_fmpq(t, v, prec + 5); + arb_mul(y, y, t, prec + 5); + arb_inv(y, y, prec); arb_neg(y, y); + arb_clear(t); fmpq_clear(v); +} + +void +arb_const_airy_bi0_eval(arb_t y, slong prec) +{ + arb_t t; fmpq_t v; arb_init(t); fmpq_init(v); + arb_set_ui(y, 3); + arb_root(y, y, 6, prec + 5); + fmpq_set_si(v, 2, 3); arb_gamma_fmpq(t, v, prec + 5); + arb_mul(y, y, t, prec + 5); + arb_inv(y, y, prec); + arb_clear(t); fmpq_clear(v); +} + +void +arb_const_airy_bi1_eval(arb_t y, slong prec) +{ + arb_t t; fmpq_t v; arb_init(t); fmpq_init(v); + arb_set_ui(y, 3); + arb_root(y, y, 6, prec + 5); + fmpq_set_si(v, 1, 3); arb_gamma_fmpq(t, v, prec + 5); + arb_div(y, y, t, prec); + arb_clear(t); fmpq_clear(v); +} + +ARB_DEF_CACHED_CONSTANT(arb_const_airy_ai0, arb_const_airy_ai0_eval) +ARB_DEF_CACHED_CONSTANT(arb_const_airy_ai1, arb_const_airy_ai1_eval) +ARB_DEF_CACHED_CONSTANT(arb_const_airy_bi0, arb_const_airy_bi0_eval) +ARB_DEF_CACHED_CONSTANT(arb_const_airy_bi1, arb_const_airy_bi1_eval) + +static void +acb_hypgeom_airy_0f1_sum_inner(acb_t s, acb_srcptr t, slong m, slong n, slong alpha, int real, slong prec) +{ + slong j, k; + mp_limb_t c, chi, clo; + + acb_zero(s); + + /* not implemented (coefficient overflow) */ + if (FLINT_BITS == 32 && n > 37000) + { + acb_indeterminate(s); + return; + } + + c = 1; + j = (n - 1) % m; + + for (k = n - 1; k >= 0; k--) + { + if (k != 0) + { + umul_ppmm(chi, clo, c, 3 * k + alpha); + + if (chi == 0) + umul_ppmm(chi, clo, clo, k); + + if (chi != 0) + { + acb_div_ui(s, s, c, prec); + c = 1; + } + } + + if (real) + arb_addmul_ui(acb_realref(s), acb_realref(t + j), c, prec); + else + acb_addmul_ui(s, t + j, c, prec); + + if (k != 0) + { + c = c * k * (3 * k + alpha); + + if (j == 0) + { + acb_mul(s, s, t + m, prec); + j = m - 1; + } + else + { + j--; + } + } + } + + acb_div_ui(s, s, c, prec); +} + +/* s1 = 0F1(1/3, z/3) + s2 = 0F1(2/3, z/3) + s4 = 0F1(4/3, z/3) + s5 = 0F1(5/3, z/3) */ +static void +acb_hypgeom_airy_0f1_sum(acb_t s1, acb_t s2, acb_t s4, acb_t s5, const acb_t z, slong n, int real, slong prec) +{ + acb_ptr t; + slong m; + + m = 2 * n_sqrt(n); + m = FLINT_MAX(m, 1); + + t = _acb_vec_init(m + 1); + _acb_vec_set_powers(t, z, m + 1, prec); + + if (s1 != NULL) acb_hypgeom_airy_0f1_sum_inner(s1, t, m, n, -2, real, prec); + if (s2 != NULL) acb_hypgeom_airy_0f1_sum_inner(s2, t, m, n, -1, real, prec); + if (s4 != NULL) acb_hypgeom_airy_0f1_sum_inner(s4, t, m, n, +1, real, prec); + if (s5 != NULL) acb_hypgeom_airy_0f1_sum_inner(s5, t, m, n, +2, real, prec); + + _acb_vec_clear(t, m + 1); +} + +void +acb_hypgeom_airy_direct(acb_t ai, acb_t aip, acb_t bi, acb_t bip, const acb_t z, slong n, slong prec) +{ + mag_t err, wmag, tm; + int is_real; + acb_t s1, s2, s4, s5, t, u; + arb_t ai0, ai1, bi0, bi1; + + mag_init(err); + mag_init(wmag); + mag_init(tm); + + acb_init(s1); + acb_init(s2); + acb_init(s4); + acb_init(s5); + acb_init(t); + acb_init(u); + + arb_init(ai0); + arb_init(ai1); + arb_init(bi0); + arb_init(bi1); + + n = FLINT_MAX(n, 2); + is_real = acb_is_real(z); + acb_get_mag(wmag, z); + + /* + With w = z^3/9, the terms are bounded by 3 w^n / [(n-1)!]^2. + + 3 w^n w w^2 + ---------- [ 1 + --- + ------- + ....] + ((n-1)!)^2 n^2 (n+1)^2 + */ + mag_pow_ui(wmag, wmag, 3); + mag_div_ui(wmag, wmag, 9); + mag_pow_ui(err, wmag, n); + + mag_div_ui(tm, wmag, n); + mag_div_ui(tm, tm, n); + mag_geom_series(tm, tm, 0); + mag_mul(err, err, tm); + + mag_rfac_ui(tm, n - 1); + mag_mul(tm, tm, tm); + mag_mul(err, err, tm); + mag_mul_ui(err, err, 3); + + acb_cube(t, z, prec); + acb_div_ui(t, t, 3, prec); + acb_hypgeom_airy_0f1_sum( + (aip != NULL || bip != NULL) ? s1 : NULL, + (ai != NULL || bi != NULL) ? s2 : NULL, + (ai != NULL || bi != NULL) ? s4 : NULL, + (aip != NULL || bip != NULL) ? s5 : NULL, t, n, is_real, prec); + + if (is_real) + { + arb_add_error_mag(acb_realref(s1), err); + arb_add_error_mag(acb_realref(s2), err); + arb_add_error_mag(acb_realref(s4), err); + arb_add_error_mag(acb_realref(s5), err); + } + else + { + acb_add_error_mag(s1, err); + acb_add_error_mag(s2, err); + acb_add_error_mag(s4, err); + acb_add_error_mag(s5, err); + } + + if (ai != NULL || aip != NULL) + { + arb_const_airy_ai0(ai0, prec); + arb_const_airy_ai1(ai1, prec); + } + + if (bi != NULL || bip != NULL) + { + arb_const_airy_bi0(bi0, prec); + arb_const_airy_bi1(bi1, prec); + } + + /* support aliasing with z */ + acb_set(t, z); + + if (ai != NULL || bi != NULL) + { + acb_mul(u, s4, t, prec); + + if (ai != NULL) + { + acb_mul_arb(ai, s2, ai0, prec); + acb_addmul_arb(ai, u, ai1, prec); + } + + if (bi != NULL) + { + acb_mul_arb(bi, s2, bi0, prec); + acb_addmul_arb(bi, u, bi1, prec); + } + } + + if (aip != NULL || bip != NULL) + { + acb_mul(u, t, t, prec); + acb_mul_2exp_si(u, u, -1); + acb_mul(u, u, s5, prec); + + if (aip != NULL) + { + acb_mul_arb(aip, s1, ai1, prec); + acb_addmul_arb(aip, u, ai0, prec); + } + + if (bip != NULL) + { + acb_mul_arb(bip, s1, bi1, prec); + acb_addmul_arb(bip, u, bi0, prec); + } + } + + mag_clear(err); + mag_clear(wmag); + mag_clear(tm); + + acb_clear(s1); + acb_clear(s2); + acb_clear(s4); + acb_clear(s5); + acb_clear(t); + acb_clear(u); + + arb_clear(ai0); + arb_clear(ai1); + arb_clear(bi0); + arb_clear(bi1); +} + diff --git a/acb_hypgeom/test/t-airy.c b/acb_hypgeom/test/t-airy.c new file mode 100644 index 00000000..7a96c225 --- /dev/null +++ b/acb_hypgeom/test/t-airy.c @@ -0,0 +1,173 @@ +/*============================================================================= + + This file is part of ARB. + + ARB is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + ARB is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ARB; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2015 Fredrik Johansson + +******************************************************************************/ + +#include "acb_hypgeom.h" + +int main() +{ + slong iter; + flint_rand_t state; + + flint_printf("airy...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 10000; iter++) + { + acb_t z, t, w; + acb_t ai1, aip1, bi1, bip1; + acb_t ai2, aip2, bi2, bip2; + slong n1, n2, prec1, prec2; + unsigned int mask; + + acb_init(z); acb_init(t); acb_init(w); + acb_init(ai1); acb_init(aip1); acb_init(bi1); acb_init(bip1); + acb_init(ai2); acb_init(aip2); acb_init(bi2); acb_init(bip2); + + prec1 = 2 + n_randint(state, 1000); + prec2 = 2 + n_randint(state, 1000); + + n1 = n_randint(state, 300); + n2 = n_randint(state, 300); + + acb_randtest_param(z, state, 1 + n_randint(state, 1000), 1 + n_randint(state, 100)); + acb_randtest_param(t, state, 1 + n_randint(state, 1000), 1 + n_randint(state, 100)); + acb_add(z, z, t, 1000); + acb_sub(z, z, t, 1000); + + switch (n_randint(state, 3)) + { + case 0: + acb_hypgeom_airy_direct(ai1, aip1, bi1, bip1, z, n1, prec1); + break; + case 1: + acb_hypgeom_airy_asymp(ai1, aip1, bi1, bip1, z, n1, prec1); + break; + default: + acb_hypgeom_airy(ai1, aip1, bi1, bip1, z, prec1); + break; + } + + switch (n_randint(state, 3)) + { + case 0: + acb_hypgeom_airy_direct(ai2, aip2, bi2, bip2, z, n2, prec2); + break; + case 1: + acb_hypgeom_airy_asymp(ai2, aip2, bi2, bip2, z, n2, prec2); + break; + default: + acb_hypgeom_airy(ai2, aip2, bi2, bip2, z, prec2); + break; + } + + if (!acb_overlaps(ai1, ai2)) + { + flint_printf("FAIL: consistency (Ai)\n\n"); + flint_printf("z = "); acb_printd(z, 30); flint_printf("\n\n"); + flint_printf("ai1 = "); acb_printd(ai1, 30); flint_printf("\n\n"); + flint_printf("ai2 = "); acb_printd(ai2, 30); flint_printf("\n\n"); + abort(); + } + + if (!acb_overlaps(aip1, aip2)) + { + flint_printf("FAIL: consistency (Ai')\n\n"); + flint_printf("z = "); acb_printd(z, 30); flint_printf("\n\n"); + flint_printf("aip1 = "); acb_printd(aip1, 30); flint_printf("\n\n"); + flint_printf("aip2 = "); acb_printd(aip2, 30); flint_printf("\n\n"); + abort(); + } + + if (!acb_overlaps(bi1, bi2)) + { + flint_printf("FAIL: consistency (Bi)\n\n"); + flint_printf("z = "); acb_printd(z, 30); flint_printf("\n\n"); + flint_printf("bi1 = "); acb_printd(bi1, 30); flint_printf("\n\n"); + flint_printf("bi2 = "); acb_printd(bi2, 30); flint_printf("\n\n"); + abort(); + } + + if (!acb_overlaps(bip1, bip2)) + { + flint_printf("FAIL: consistency (Bi')\n\n"); + flint_printf("z = "); acb_printd(z, 30); flint_printf("\n\n"); + flint_printf("bip1 = "); acb_printd(bip1, 30); flint_printf("\n\n"); + flint_printf("bip2 = "); acb_printd(bip2, 30); flint_printf("\n\n"); + abort(); + } + + acb_mul(w, ai1, bip1, prec1); + acb_submul(w, bi1, aip1, prec1); + acb_const_pi(t, prec1); + acb_inv(t, t, prec1); + + if (!acb_overlaps(w, t)) + { + flint_printf("FAIL: wronskian\n\n"); + flint_printf("z = "); acb_printd(z, 30); flint_printf("\n\n"); + flint_printf("ai1 = "); acb_printd(ai1, 30); flint_printf("\n\n"); + flint_printf("aip1 = "); acb_printd(aip1, 30); flint_printf("\n\n"); + flint_printf("bi1 = "); acb_printd(bi1, 30); flint_printf("\n\n"); + flint_printf("bip1 = "); acb_printd(bip1, 30); flint_printf("\n\n"); + flint_printf("w = "); acb_printd(w, 30); flint_printf("\n\n"); + abort(); + } + + mask = n_randlimb(state); + + acb_hypgeom_airy((mask & 1) ? ai2 : NULL, + (mask & 2) ? aip2 : NULL, + (mask & 4) ? bi2 : NULL, + (mask & 8) ? bip2 : NULL, z, prec2); + + if (!acb_overlaps(ai1, ai2) || !acb_overlaps(aip1, aip2) || + !acb_overlaps(bi1, bi2) || !acb_overlaps(bip1, bip2)) + { + flint_printf("FAIL: consistency (mask)\n\n"); + flint_printf("mask = %u\n\n", mask); + flint_printf("z = "); acb_printd(z, 30); flint_printf("\n\n"); + flint_printf("ai1 = "); acb_printd(ai1, 30); flint_printf("\n\n"); + flint_printf("ai2 = "); acb_printd(ai2, 30); flint_printf("\n\n"); + flint_printf("aip1 = "); acb_printd(aip1, 30); flint_printf("\n\n"); + flint_printf("aip2 = "); acb_printd(aip2, 30); flint_printf("\n\n"); + flint_printf("bi1 = "); acb_printd(bi1, 30); flint_printf("\n\n"); + flint_printf("bi2 = "); acb_printd(bi2, 30); flint_printf("\n\n"); + flint_printf("bip1 = "); acb_printd(bip1, 30); flint_printf("\n\n"); + flint_printf("bip2 = "); acb_printd(bip2, 30); flint_printf("\n\n"); + abort(); + } + + acb_clear(z); acb_clear(t); acb_clear(w); + acb_clear(ai1); acb_clear(aip1); acb_clear(bi1); acb_clear(bip1); + acb_clear(ai2); acb_clear(aip2); acb_clear(bi2); acb_clear(bip2); + } + + flint_randclear(state); + flint_cleanup(); + flint_printf("PASS\n"); + return EXIT_SUCCESS; +} diff --git a/acb_hypgeom/test/t-airy_bound.c b/acb_hypgeom/test/t-airy_bound.c new file mode 100644 index 00000000..bcb4c7b3 --- /dev/null +++ b/acb_hypgeom/test/t-airy_bound.c @@ -0,0 +1,158 @@ +/*============================================================================= + + This file is part of ARB. + + ARB is free software; you can redistribute it and/or modify + it under the terms of the GNU General Public License as published by + the Free Software Foundation; either version 2 of the License, or + (at your option) any later version. + + ARB is distributed in the hope that it will be useful, + but WITHOUT ANY WARRANTY; without even the implied warranty of + MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + GNU General Public License for more details. + + You should have received a copy of the GNU General Public License + along with ARB; if not, write to the Free Software + Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA + +=============================================================================*/ +/****************************************************************************** + + Copyright (C) 2015 Fredrik Johansson + +******************************************************************************/ + +#include "acb_hypgeom.h" + +int main() +{ + slong iter; + flint_rand_t state; + + flint_printf("airy_bound...."); + fflush(stdout); + + flint_randinit(state); + + for (iter = 0; iter < 10000; iter++) + { + acb_t ai, aip, bi, bip, z1, z2; + slong prec; + mag_t aib, aipb, bib, bipb, aim, aipm, bim, bipm; + + acb_init(ai); + acb_init(aip); + acb_init(bi); + acb_init(bip); + acb_init(z1); + acb_init(z2); + + mag_init(aib); + mag_init(aipb); + mag_init(bib); + mag_init(bipb); + + mag_init(aim); + mag_init(aipm); + mag_init(bim); + mag_init(bipm); + + acb_randtest(z1, state, 1 + n_randint(state, 200), 1 + n_randint(state, 10)); + arb_mul_ui(acb_realref(z1), acb_realref(z1), n_randint(state, 300), 1 + n_randint(state, 200)); + arb_mul_ui(acb_imagref(z1), acb_imagref(z1), n_randint(state, 300), 1 + n_randint(state, 200)); + + acb_zero(z2); + + arf_set_mag(arb_midref(acb_realref(z2)), arb_radref(acb_realref(z1))); + arf_set_mag(arb_midref(acb_imagref(z2)), arb_radref(acb_imagref(z1))); + + switch (n_randint(state, 5)) + { + case 0: + arf_add(arb_midref(acb_realref(z2)), arb_midref(acb_realref(z1)), arb_midref(acb_realref(z2)), ARF_PREC_EXACT, ARF_RND_DOWN); + arf_add(arb_midref(acb_imagref(z2)), arb_midref(acb_imagref(z1)), arb_midref(acb_imagref(z2)), ARF_PREC_EXACT, ARF_RND_DOWN); + break; + case 1: + arf_add(arb_midref(acb_realref(z2)), arb_midref(acb_realref(z1)), arb_midref(acb_realref(z2)), ARF_PREC_EXACT, ARF_RND_DOWN); + arf_sub(arb_midref(acb_imagref(z2)), arb_midref(acb_imagref(z1)), arb_midref(acb_imagref(z2)), ARF_PREC_EXACT, ARF_RND_DOWN); + break; + case 2: + arf_sub(arb_midref(acb_realref(z2)), arb_midref(acb_realref(z1)), arb_midref(acb_realref(z2)), ARF_PREC_EXACT, ARF_RND_DOWN); + arf_add(arb_midref(acb_imagref(z2)), arb_midref(acb_imagref(z1)), arb_midref(acb_imagref(z2)), ARF_PREC_EXACT, ARF_RND_DOWN); + break; + case 3: + arf_sub(arb_midref(acb_realref(z2)), arb_midref(acb_realref(z1)), arb_midref(acb_realref(z2)), ARF_PREC_EXACT, ARF_RND_DOWN); + arf_sub(arb_midref(acb_imagref(z2)), arb_midref(acb_imagref(z1)), arb_midref(acb_imagref(z2)), ARF_PREC_EXACT, ARF_RND_DOWN); + break; + default: + arf_set(arb_midref(acb_realref(z2)), arb_midref(acb_realref(z1))); + arf_set(arb_midref(acb_imagref(z2)), arb_midref(acb_imagref(z1))); + } + + acb_hypgeom_airy_bound(aib, aipb, bib, bipb, z1); + + prec = MAG_BITS + 10; + + do { + acb_hypgeom_airy(ai, aip, bi, bip, z2, prec); + + if (acb_rel_accuracy_bits(ai) >= MAG_BITS && + acb_rel_accuracy_bits(aip) >= MAG_BITS && + acb_rel_accuracy_bits(bi) >= MAG_BITS && + acb_rel_accuracy_bits(bip) >= MAG_BITS) + break; + + prec *= 2; + } while (1); + + acb_get_mag(aim, ai); + acb_get_mag(aipm, aip); + acb_get_mag(bim, bi); + acb_get_mag(bipm, bip); + + if (mag_cmp(aim, aib) > 0 || mag_cmp(aipm, aipb) > 0 || + mag_cmp(bim, bib) > 0 || mag_cmp(aipm, bipb) > 0) + { + printf("FAIL\n"); + flint_printf("z1 = "); acb_printd(z1, 20); flint_printf("\n"); + flint_printf("z2 = "); acb_printd(z2, 20); flint_printf("\n\n"); + flint_printf("ai = "); acb_printd(ai, 20); flint_printf("\n"); + flint_printf("aim = "); mag_printd(aim, 10); printf("\n"); + flint_printf("aib = "); mag_printd(aib, 10); printf("\n\n"); + flint_printf("api = "); acb_printd(aip, 20); flint_printf("\n"); + flint_printf("aipm = "); mag_printd(aipm, 10); printf("\n"); + flint_printf("aipb = "); mag_printd(aipb, 10); printf("\n\n"); + flint_printf("bi = "); acb_printd(bi, 20); flint_printf("\n"); + flint_printf("bim = "); mag_printd(bim, 10); printf("\n"); + flint_printf("bib = "); mag_printd(bib, 10); printf("\n\n"); + flint_printf("bpi = "); acb_printd(bip, 20); flint_printf("\n"); + flint_printf("bipm = "); mag_printd(bipm, 10); printf("\n"); + flint_printf("bipb = "); mag_printd(bipb, 10); printf("\n\n"); + abort(); + } + + acb_clear(ai); + acb_clear(aip); + acb_clear(bi); + acb_clear(bip); + acb_clear(z1); + acb_clear(z2); + + mag_clear(aib); + mag_clear(aipb); + mag_clear(bib); + mag_clear(bipb); + + mag_clear(aim); + mag_clear(aipm); + mag_clear(bim); + mag_clear(bipm); + } + + flint_randclear(state); + flint_cleanup(); + flint_printf("PASS\n"); + return EXIT_SUCCESS; +} + diff --git a/doc/source/acb_hypgeom.rst b/doc/source/acb_hypgeom.rst index fe9cbf3c..1297a5d0 100644 --- a/doc/source/acb_hypgeom.rst +++ b/doc/source/acb_hypgeom.rst @@ -475,6 +475,148 @@ Bessel functions Computes the modified Bessel function of the second kind `K_{\nu}(z)` using an automatic algorithm choice. +Airy functions +------------------------------------------------------------------------------- + +The Airy functions are linearly independent solutions of the +differential equation `y'' - zy = 0`. All solutions are entire functions. +The standard solutions are denoted `\operatorname{Ai}(z), \operatorname{Bi}(z)`. +For negative *z*, both functions are oscillatory. For positive *z*, +the first function decreases exponentially while the second increases +exponentially. + +The Airy functions can be expressed in terms of Bessel functions of fractional +order, but this is inconvenient since such formulas +only hold piecewise (due to the Stokes phenomenon). Computation of the +Airy functions can also be optimized more than Bessel functions in general. +We therefore provide a dedicated interface for evaluating Airy functions. + +The following functions optionally compute +`(\operatorname{Ai}(z), \operatorname{Ai}'(z), \operatorname{Bi}(z), \operatorname{Bi}'(z))` +simultaneously. Any of the four functions can be omitted by passing +*NULL* for the unwanted output variables. +Note that higher derivatives of the Airy functions can be computed +via recurrence relations. + +.. function:: void acb_hypgeom_airy_direct(acb_t ai, acb_t ai_prime, acb_t bi, acb_t bi_prime, const acb_t z, slong n, slong prec) + + Computes the Airy functions using direct series expansions truncated at *n* terms. + Error bounds are included in the output. + +.. function:: void acb_hypgeom_airy_asymp(acb_t ai, acb_t ai_prime, acb_t bi, acb_t bi_prime, const acb_t z, slong n, slong prec) + + Computes the Airy functions using asymptotic expansions truncated at *n* terms. + Error bounds, based on Olver (DLMF section 9.7), are included in the output. + For `\arg(z) < \pi` and `\zeta = (2/3) z^{3/2}`, we have + + .. math :: + + \operatorname{Ai}(z) = \frac{e^{-\zeta}}{2 \sqrt{\pi} z^{1/4}} \left[S_n(\zeta) + R_n(\zeta)\right], \quad + \operatorname{Ai}'(z) = -\frac{z^{1/4} e^{-\zeta}}{2 \sqrt{\pi}} \left[(S'_n(\zeta) + R'_n(\zeta)\right] + + S_n(\zeta) = \sum_{k=0}^{n-1} (-1)^k \frac{u(k)}{\zeta^k}, \quad + S'_n(\zeta) = \sum_{k=0}^{n-1} (-1)^k \frac{v(k)}{\zeta^k} + + u(k) = \frac{(1/6)_k (5/6)_k}{2^k k!}, \quad + v(k) = \frac{6k+1}{1-6k} u(k). + + Assuming that *n* is positive, the error terms are bounded by + + .. math :: + + |R_n(\zeta)| \le C |u(n)| |\zeta|^{-n}, \quad |R'_n(\zeta)| \le C |v(n)| |\zeta|^{-n} + + where + + .. math :: + + C = \begin{cases} + 2 \exp(7 / (36 |\zeta|)) & |\arg(z)| \le \pi/3 \\ + 2 \chi(n) \exp(7 \pi / (72 |\zeta|)) & \pi/3 \le |\arg(z)| \le 2\pi/3 \\ + 4 \chi(n) \exp(7 \pi / (36 |\operatorname{re}(\zeta)|)) |\cos(\arg(\zeta))|^{-n} & 2\pi/3 \le |\arg(z)| < \pi. + \end{cases} + + For computing Bi when *z* is roughly in the positive half-plane, we use the + connection formulas + + .. math :: + + \operatorname{Bi}(z) = -i (2 w^{+1} \operatorname{Ai}(z w^{-2}) - \operatorname{Ai}(z)) + + \operatorname{Bi}(z) = +i (2 w^{-1} \operatorname{Ai}(z w^{+2}) - \operatorname{Ai}(z)) + + where `w = \exp(\pi i/3)`. Combining roots of unity gives + + .. math :: + + \operatorname{Bi}(z) = \frac{1}{2 \sqrt{\pi} z^{1/4}} [2X + iY] + + \operatorname{Bi}(z) = \frac{1}{2 \sqrt{\pi} z^{1/4}} [2X - iY] + + X = \exp(+\zeta) [S_n(-\zeta) + R_n(-\zeta)], \quad Y = \exp(-\zeta) [S_n(\zeta) + R_n(\zeta)] + + where the upper formula is valid + for `-\pi/3 < \arg(z) < \pi` and the lower formula is valid for `-\pi < \arg(z) < \pi/3`. + We proceed analogously for the derivative of Bi. + + In the negative half-plane, we use the connection formulas + + .. math :: + + \operatorname{Ai}(z) = e^{+\pi i/3} \operatorname{Ai}(z_1) + e^{-\pi i/3} \operatorname{Ai}(z_2) + + \operatorname{Bi}(z) = e^{-\pi i/6} \operatorname{Ai}(z_1) + e^{+\pi i/6} \operatorname{Ai}(z_2) + + where `z_1 = -z e^{+\pi i/3}`, `z_2 = -z e^{-\pi i/3}`. + Provided that `|\arg(-z)| < 2 \pi / 3`, we have + `|\arg(z_1)|, |\arg(z_2)| < \pi`, and thus the asymptotic expansion + for Ai can be used. As before, we collect roots of unity to obtain + + .. math :: + + \operatorname{Ai}(z) = A_1 [S_n(i \zeta) + R_n(i \zeta)] + + A_2 [S_n(-i \zeta) + R_n(-i \zeta)] + + \operatorname{Bi}(z) = A_3 [S_n(i \zeta) + R_n(i \zeta)] + + A_4 [S_n(-i \zeta) + R_n(-i \zeta)] + + where `\zeta = (2/3) (-z)^{3/2}` and + + .. math :: + + A_1 = \frac{\exp(-i (\zeta - \pi/4))}{2 \sqrt{\pi} (-z)^{1/4}}, \quad + A_2 = \frac{\exp(+i (\zeta - \pi/4))}{2 \sqrt{\pi} (-z)^{1/4}}, \quad + A_3 = -i A_1, \quad + A_4 = +i A_2. + + The differentiated formulas are analogous. + +.. function:: void acb_hypgeom_airy_bound(mag_t ai, mag_t aip, mag_t bi, mag_t bip, const acb_t z) + + Computes bounds for the Airy functions using first-order asymptotic + expansions together with error bounds. This function uses some + shortcuts to make it slightly faster than calling + :func:`acb_hypgeom_airy_asymp` with `n = 1`. + +.. function:: void acb_hypgeom_airy(acb_t ai, acb_t ai_prime, acb_t bi, acb_t bi_prime, const acb_t z, slong prec) + + Computes Airy functions using an automatic algorithm choice. + + We use :func:`acb_hypgeom_airy_asymp` whenever this gives full accuracy + and :func:`acb_hypgeom_airy_direct` otherwise. + + In the latter case, we first use double precision arithmetic to + determine an accurate estimate of the working precision needed + to compute the Airy functions accurately for given *z*. This estimate is + obtained by comparing the leading-order asymptotic estimate of the Airy + functions with the magnitude of the largest term in the power series. + The estimate is generic in the sense that it does not take into account + vanishing near the roots of the functions. + + We subsequently evaluate the power series at the midpoint of *z* and + bound the propagated error using derivatives. Derivatives are + bounded using :func:`acb_hypgeom_airy_bound`. + Incomplete gamma functions -------------------------------------------------------------------------------