hiro/arb
1
0
Fork 0
mirror of https://github.com/vale981/arb synced 2025-03-05 09:21:38 -05:00
Code Issues Projects Releases Packages Wiki Activity Actions

initial code for Airy functions

Browse source
This commit is contained in:
Fredrik Johansson 2015-11-19 11:10:22 +01:00
parent f3beae93aa
commit 9eed54684c
8 changed files with 1926 additions and 0 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

5
acb_hypgeom.h
View file

@ -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);

294
acb_hypgeom/airy.c Normal file
View file

@ -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);
}

601
acb_hypgeom/airy_asymp.c Normal file
View file

@ -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);
}

252
acb_hypgeom/airy_bound.c Normal file
View file

@ -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);
}

301
acb_hypgeom/airy_direct.c Normal file
View file

@ -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);
}

173
acb_hypgeom/test/t-airy.c Normal file
View file

@ -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;
}

158
acb_hypgeom/test/t-airy_bound.c Normal file
View file

@ -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;
}

142
doc/source/acb_hypgeom.rst
View file

@ -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
-------------------------------------------------------------------------------