From 7f98274601618bc2c0d783b62f4a582f2cd6a3eb Mon Sep 17 00:00:00 2001 From: fredrik Date: Fri, 23 Mar 2018 13:41:46 +0100 Subject: [PATCH] implement scaled modified Bessel functions --- acb_hypgeom.h | 12 ++++---- acb_hypgeom/0f1.c | 2 +- acb_hypgeom/bessel_i.c | 56 +++++++++++++++++++++++++++++------ acb_hypgeom/bessel_k.c | 56 ++++++++++++++++++++++++++++------- acb_hypgeom/test/t-bessel_i.c | 54 +++++++++++++++++++++++---------- acb_hypgeom/test/t-bessel_k.c | 48 ++++++++++++++++++++++-------- arb_hypgeom.h | 3 ++ arb_hypgeom/test/t-wrappers.c | 8 +++++ arb_hypgeom/wrappers.c | 34 +++++++++++++++++++++ doc/source/acb_hypgeom.rst | 28 ++++++++++++++---- doc/source/arb_hypgeom.rst | 8 +++++ 11 files changed, 250 insertions(+), 59 deletions(-) diff --git a/acb_hypgeom.h b/acb_hypgeom.h index 591ff256..f0bded86 100644 --- a/acb_hypgeom.h +++ b/acb_hypgeom.h @@ -111,14 +111,16 @@ void acb_hypgeom_bessel_j_0f1(acb_t res, const acb_t nu, const acb_t z, slong pr void acb_hypgeom_bessel_j_asymp(acb_t res, const acb_t nu, const acb_t z, slong prec); void acb_hypgeom_bessel_j(acb_t res, const acb_t nu, const acb_t z, slong prec); -void acb_hypgeom_bessel_i_0f1(acb_t res, const acb_t nu, const acb_t z, slong prec); -void acb_hypgeom_bessel_i_asymp(acb_t res, const acb_t nu, const acb_t z, slong prec); +void acb_hypgeom_bessel_i_0f1(acb_t res, const acb_t nu, const acb_t z, int scaled, slong prec); +void acb_hypgeom_bessel_i_asymp(acb_t res, const acb_t nu, const acb_t z, int scaled, slong prec); void acb_hypgeom_bessel_i(acb_t res, const acb_t nu, const acb_t z, slong prec); +void acb_hypgeom_bessel_i_scaled(acb_t res, const acb_t nu, const acb_t z, slong prec); -void acb_hypgeom_bessel_k_0f1(acb_t res, const acb_t nu, const acb_t z, slong prec); -void acb_hypgeom_bessel_k_0f1_series(acb_poly_t res, const acb_poly_t n, const acb_poly_t z, slong len, slong prec); -void acb_hypgeom_bessel_k_asymp(acb_t res, const acb_t nu, const acb_t z, slong prec); +void acb_hypgeom_bessel_k_0f1(acb_t res, const acb_t nu, const acb_t z, int scaled, slong prec); +void acb_hypgeom_bessel_k_0f1_series(acb_poly_t res, const acb_poly_t n, const acb_poly_t z, int scaled, slong len, slong prec); +void acb_hypgeom_bessel_k_asymp(acb_t res, const acb_t nu, const acb_t z, int scaled, slong prec); void acb_hypgeom_bessel_k(acb_t res, const acb_t nu, const acb_t z, slong prec); +void acb_hypgeom_bessel_k_scaled(acb_t res, const acb_t nu, const acb_t z, slong prec); void acb_hypgeom_bessel_y(acb_t res, const acb_t nu, const acb_t z, slong prec); void acb_hypgeom_bessel_jy(acb_t res1, acb_t res2, const acb_t nu, const acb_t z, slong prec); diff --git a/acb_hypgeom/0f1.c b/acb_hypgeom/0f1.c index e2c063b8..5134d4c2 100644 --- a/acb_hypgeom/0f1.c +++ b/acb_hypgeom/0f1.c @@ -37,7 +37,7 @@ acb_hypgeom_0f1_asymp(acb_t res, const acb_t a, const acb_t z, int regularized, if (neg) acb_hypgeom_bessel_j_asymp(v, u, v, prec); else - acb_hypgeom_bessel_i_asymp(v, u, v, prec); + acb_hypgeom_bessel_i_asymp(v, u, v, 0, prec); acb_neg(u, u); acb_pow(t, t, u, prec); diff --git a/acb_hypgeom/bessel_i.c b/acb_hypgeom/bessel_i.c index 898cfe23..239f8ef7 100644 --- a/acb_hypgeom/bessel_i.c +++ b/acb_hypgeom/bessel_i.c @@ -13,7 +13,7 @@ void acb_hypgeom_bessel_i_asymp_prefactors(acb_t A, acb_t B, acb_t C, - const acb_t nu, const acb_t z, slong prec) + const acb_t nu, const acb_t z, int scaled, slong prec) { acb_t t, u; @@ -53,15 +53,26 @@ acb_hypgeom_bessel_i_asymp_prefactors(acb_t A, acb_t B, acb_t C, arb_union(acb_imagref(t), acb_imagref(t), acb_imagref(u), prec); } - acb_exp_invexp(B, A, z, prec); - acb_mul(A, A, t, prec); + if (scaled) + { + acb_neg(u, z); + acb_mul_2exp_si(u, u, 1); + acb_exp(u, u, prec); + acb_mul(A, t, u, prec); + acb_one(B); + } + else + { + acb_exp_invexp(B, A, z, prec); + acb_mul(A, A, t, prec); + } acb_clear(t); acb_clear(u); } void -acb_hypgeom_bessel_i_asymp(acb_t res, const acb_t nu, const acb_t z, slong prec) +acb_hypgeom_bessel_i_asymp(acb_t res, const acb_t nu, const acb_t z, int scaled, slong prec) { acb_t A1, A2, C, U1, U2, s, t, u; int is_real, is_imag; @@ -89,7 +100,10 @@ acb_hypgeom_bessel_i_asymp(acb_t res, const acb_t nu, const acb_t z, slong prec) is_imag = 1; } - acb_hypgeom_bessel_i_asymp_prefactors(A1, A2, C, nu, z, prec); + if (scaled) + is_imag = 0; + + acb_hypgeom_bessel_i_asymp_prefactors(A1, A2, C, nu, z, scaled, prec); /* todo: if Ap ~ 2^a and Am = 2^b and U1 ~ U2 ~ 1, change precision? */ @@ -134,7 +148,7 @@ acb_hypgeom_bessel_i_asymp(acb_t res, const acb_t nu, const acb_t z, slong prec) } void -acb_hypgeom_bessel_i_0f1(acb_t res, const acb_t nu, const acb_t z, slong prec) +acb_hypgeom_bessel_i_0f1(acb_t res, const acb_t nu, const acb_t z, int scaled, slong prec) { acb_struct b[2]; acb_t w, c, t; @@ -143,7 +157,7 @@ acb_hypgeom_bessel_i_0f1(acb_t res, const acb_t nu, const acb_t z, slong prec) { acb_init(t); acb_neg(t, nu); - acb_hypgeom_bessel_i_0f1(res, t, z, prec); + acb_hypgeom_bessel_i_0f1(res, t, z, scaled, prec); acb_clear(t); return; } @@ -169,6 +183,13 @@ acb_hypgeom_bessel_i_0f1(acb_t res, const acb_t nu, const acb_t z, slong prec) acb_hypgeom_pfq_direct(t, NULL, 0, b, 2, w, -1, prec); + if (scaled) + { + acb_neg(w, z); + acb_exp(w, w, prec); + acb_mul(t, t, w, prec); + } + acb_mul(res, t, c, prec); acb_clear(b + 0); @@ -188,9 +209,26 @@ acb_hypgeom_bessel_i(acb_t res, const acb_t nu, const acb_t z, slong prec) if (mag_cmp_2exp_si(zmag, 4) < 0 || (mag_cmp_2exp_si(zmag, 64) < 0 && 2 * mag_get_d(zmag) < prec)) - acb_hypgeom_bessel_i_0f1(res, nu, z, prec); + acb_hypgeom_bessel_i_0f1(res, nu, z, 0, prec); else - acb_hypgeom_bessel_i_asymp(res, nu, z, prec); + acb_hypgeom_bessel_i_asymp(res, nu, z, 0, prec); + + mag_clear(zmag); +} + +void +acb_hypgeom_bessel_i_scaled(acb_t res, const acb_t nu, const acb_t z, slong prec) +{ + mag_t zmag; + + mag_init(zmag); + acb_get_mag(zmag, z); + + if (mag_cmp_2exp_si(zmag, 4) < 0 || + (mag_cmp_2exp_si(zmag, 64) < 0 && 2 * mag_get_d(zmag) < prec)) + acb_hypgeom_bessel_i_0f1(res, nu, z, 1, prec); + else + acb_hypgeom_bessel_i_asymp(res, nu, z, 1, prec); mag_clear(zmag); } diff --git a/acb_hypgeom/bessel_k.c b/acb_hypgeom/bessel_k.c index c8cd0196..83f7792c 100644 --- a/acb_hypgeom/bessel_k.c +++ b/acb_hypgeom/bessel_k.c @@ -12,7 +12,7 @@ #include "acb_hypgeom.h" void -acb_hypgeom_bessel_k_asymp(acb_t res, const acb_t nu, const acb_t z, slong prec) +acb_hypgeom_bessel_k_asymp(acb_t res, const acb_t nu, const acb_t z, int scaled, slong prec) { acb_t t, a, b, w; @@ -32,9 +32,12 @@ acb_hypgeom_bessel_k_asymp(acb_t res, const acb_t nu, const acb_t z, slong prec) acb_hypgeom_u_asymp(t, a, b, w, -1, prec); - acb_neg(w, z); - acb_exp(w, w, prec); - acb_mul(t, t, w, prec); + if (!scaled) + { + acb_neg(a, z); + acb_exp(a, a, prec); + acb_mul(t, t, a, prec); + } acb_mul_2exp_si(w, z, 1); acb_rsqrt(w, w, prec); @@ -52,7 +55,7 @@ acb_hypgeom_bessel_k_asymp(acb_t res, const acb_t nu, const acb_t z, slong prec) void acb_hypgeom_bessel_k_0f1_series(acb_poly_t res, const acb_poly_t nu, const acb_poly_t z, - slong len, slong prec) + int scaled, slong len, slong prec) { acb_poly_t s, u, A, B; acb_poly_struct b[2]; @@ -101,6 +104,12 @@ acb_hypgeom_bessel_k_0f1_series(acb_poly_t res, acb_poly_shift_right(B, B, 1); } + if (scaled) + { + acb_poly_exp_series(u, z, len, prec); + acb_poly_mullow(A, A, u, len, prec); + } + acb_poly_div_series(res, A, B, len, prec); arb_const_pi(c, prec); @@ -116,7 +125,7 @@ acb_hypgeom_bessel_k_0f1_series(acb_poly_t res, } void -acb_hypgeom_bessel_k_0f1(acb_t res, const acb_t nu, const acb_t z, slong prec) +acb_hypgeom_bessel_k_0f1(acb_t res, const acb_t nu, const acb_t z, int scaled, slong prec) { if (acb_is_int(nu)) { @@ -130,7 +139,7 @@ acb_hypgeom_bessel_k_0f1(acb_t res, const acb_t nu, const acb_t z, slong prec) acb_poly_set_coeff_si(nux, 1, 1); acb_poly_set_acb(zx, z); - acb_hypgeom_bessel_k_0f1_series(rx, nux, zx, 1, prec); + acb_hypgeom_bessel_k_0f1_series(rx, nux, zx, scaled, 1, prec); acb_poly_get_coeff_acb(res, rx, 0); @@ -176,8 +185,16 @@ acb_hypgeom_bessel_k_0f1(acb_t res, const acb_t nu, const acb_t z, slong prec) acb_mul(t, t, nu, prec); acb_div(u, u, t, prec); - acb_sub(res, v, u, prec); - acb_mul_2exp_si(res, res, -1); + acb_sub(v, v, u, prec); + acb_mul_2exp_si(v, v, -1); + + if (scaled) + { + acb_exp(t, z, prec); + acb_mul(v, v, t, prec); + } + + acb_set(res, v); acb_clear(t); acb_clear(u); @@ -198,9 +215,26 @@ acb_hypgeom_bessel_k(acb_t res, const acb_t nu, const acb_t z, slong prec) if (mag_cmp_2exp_si(zmag, 4) < 0 || (mag_cmp_2exp_si(zmag, 64) < 0 && 2 * mag_get_d(zmag) < prec)) - acb_hypgeom_bessel_k_0f1(res, nu, z, prec); + acb_hypgeom_bessel_k_0f1(res, nu, z, 0, prec); else - acb_hypgeom_bessel_k_asymp(res, nu, z, prec); + acb_hypgeom_bessel_k_asymp(res, nu, z, 0, prec); + + mag_clear(zmag); +} + +void +acb_hypgeom_bessel_k_scaled(acb_t res, const acb_t nu, const acb_t z, slong prec) +{ + mag_t zmag; + + mag_init(zmag); + acb_get_mag(zmag, z); + + if (mag_cmp_2exp_si(zmag, 4) < 0 || + (mag_cmp_2exp_si(zmag, 64) < 0 && 2 * mag_get_d(zmag) < prec)) + acb_hypgeom_bessel_k_0f1(res, nu, z, 1, prec); + else + acb_hypgeom_bessel_k_asymp(res, nu, z, 1, prec); mag_clear(zmag); } diff --git a/acb_hypgeom/test/t-bessel_i.c b/acb_hypgeom/test/t-bessel_i.c index e4df4533..af4a4388 100644 --- a/acb_hypgeom/test/t-bessel_i.c +++ b/acb_hypgeom/test/t-bessel_i.c @@ -25,6 +25,7 @@ int main() { acb_t nu, z, jv, iv, t; slong prec; + int scaled; acb_init(nu); acb_init(z); @@ -33,6 +34,7 @@ int main() acb_init(t); prec = 2 + n_randint(state, 500); + scaled = n_randint(state, 2); acb_randtest_param(nu, state, 1 + n_randint(state, 1000), 1 + n_randint(state, 10)); acb_randtest(z, state, 1 + n_randint(state, 1000), 1 + n_randint(state, 100)); @@ -40,13 +42,16 @@ int main() switch (n_randint(state, 3)) { case 0: - acb_hypgeom_bessel_i_asymp(iv, nu, z, prec); + acb_hypgeom_bessel_i_asymp(iv, nu, z, scaled, prec); break; case 1: - acb_hypgeom_bessel_i_0f1(iv, nu, z, prec); + acb_hypgeom_bessel_i_0f1(iv, nu, z, scaled, prec); break; default: - acb_hypgeom_bessel_i(iv, nu, z, prec); + if (scaled) + acb_hypgeom_bessel_i_scaled(iv, nu, z, prec); + else + acb_hypgeom_bessel_i(iv, nu, z, prec); } acb_mul_onei(t, z); @@ -57,9 +62,17 @@ int main() acb_pow(t, t, nu, prec); acb_div(jv, jv, t, prec); + if (scaled) + { + acb_neg(t, z); + acb_exp(t, t, prec); + acb_mul(jv, jv, t, prec); + } + if (!acb_overlaps(iv, jv)) { flint_printf("FAIL: consistency with bessel_j\n\n"); + flint_printf("scaled = %d\n\n", scaled); flint_printf("nu = "); acb_printd(nu, 30); flint_printf("\n\n"); flint_printf("z = "); acb_printd(z, 30); flint_printf("\n\n"); flint_printf("iv = "); acb_printd(iv, 30); flint_printf("\n\n"); @@ -105,10 +118,10 @@ int main() switch (n_randint(state, 3)) { case 0: - acb_hypgeom_bessel_i_asymp(w0, nu0, z, prec0); + acb_hypgeom_bessel_i_asymp(w0, nu0, z, 0, prec0); break; case 1: - acb_hypgeom_bessel_i_0f1(w0, nu0, z, prec0); + acb_hypgeom_bessel_i_0f1(w0, nu0, z, 0, prec0); break; default: acb_hypgeom_bessel_i(w0, nu0, z, prec0); @@ -117,10 +130,10 @@ int main() switch (n_randint(state, 3)) { case 0: - acb_hypgeom_bessel_i_asymp(w1, nu0, z, prec1); + acb_hypgeom_bessel_i_asymp(w1, nu0, z, 0, prec1); break; case 1: - acb_hypgeom_bessel_i_0f1(w1, nu0, z, prec1); + acb_hypgeom_bessel_i_0f1(w1, nu0, z, 0, prec1); break; default: acb_hypgeom_bessel_i(w1, nu0, z, prec1); @@ -139,10 +152,10 @@ int main() switch (n_randint(state, 3)) { case 0: - acb_hypgeom_bessel_i_asymp(w1, nu1, z, prec1); + acb_hypgeom_bessel_i_asymp(w1, nu1, z, 0, prec1); break; case 1: - acb_hypgeom_bessel_i_0f1(w1, nu1, z, prec1); + acb_hypgeom_bessel_i_0f1(w1, nu1, z, 0, prec1); break; default: acb_hypgeom_bessel_i(w1, nu1, z, prec1); @@ -151,13 +164,22 @@ int main() switch (n_randint(state, 3)) { case 0: - acb_hypgeom_bessel_i_asymp(w2, nu2, z, prec2); + acb_hypgeom_bessel_i_asymp(w2, nu2, z, 0, prec2); break; case 1: - acb_hypgeom_bessel_i_0f1(w2, nu2, z, prec2); + acb_hypgeom_bessel_i_0f1(w2, nu2, z, 0, prec2); break; default: - acb_hypgeom_bessel_i(w2, nu2, z, prec2); + if (n_randint(state, 2)) + { + acb_hypgeom_bessel_i(w2, nu2, z, prec2); + } + else + { + acb_hypgeom_bessel_i_scaled(w2, nu2, z, prec2); + acb_exp(t, z, prec2); + acb_mul(w2, w2, t, prec2); + } } acb_mul(t, w1, nu1, prec0); @@ -182,10 +204,10 @@ int main() switch (n_randint(state, 3)) { case 0: - acb_hypgeom_bessel_i_asymp(w2, t, z, prec2); + acb_hypgeom_bessel_i_asymp(w2, t, z, 0, prec2); break; case 1: - acb_hypgeom_bessel_i_0f1(w2, t, z, prec2); + acb_hypgeom_bessel_i_0f1(w2, t, z, 0, prec2); break; default: acb_hypgeom_bessel_i(w2, t, z, prec2); @@ -197,10 +219,10 @@ int main() switch (n_randint(state, 3)) { case 0: - acb_hypgeom_bessel_i_asymp(w2, t, z, prec2); + acb_hypgeom_bessel_i_asymp(w2, t, z, 0, prec2); break; case 1: - acb_hypgeom_bessel_i_0f1(w2, t, z, prec2); + acb_hypgeom_bessel_i_0f1(w2, t, z, 0, prec2); break; default: acb_hypgeom_bessel_i(w2, t, z, prec2); diff --git a/acb_hypgeom/test/t-bessel_k.c b/acb_hypgeom/test/t-bessel_k.c index 73d08821..91916012 100644 --- a/acb_hypgeom/test/t-bessel_k.c +++ b/acb_hypgeom/test/t-bessel_k.c @@ -25,6 +25,7 @@ int main() { acb_t nu0, nu1, nu2, z, w0, w1, w2, t, u; slong prec0, prec1, prec2; + int scaled; acb_init(nu0); acb_init(nu1); @@ -40,6 +41,7 @@ int main() prec1 = 2 + n_randint(state, 700); prec2 = 2 + n_randint(state, 700); + scaled = n_randint(state, 2); acb_randtest_param(nu0, state, 1 + n_randint(state, 1000), 1 + n_randint(state, 100)); acb_randtest(z, state, 1 + n_randint(state, 1000), 1 + n_randint(state, 100)); acb_randtest(w0, state, 1 + n_randint(state, 1000), 1 + n_randint(state, 100)); @@ -55,30 +57,37 @@ int main() switch (n_randint(state, 3)) { case 0: - acb_hypgeom_bessel_k_asymp(w0, nu0, z, prec0); + acb_hypgeom_bessel_k_asymp(w0, nu0, z, scaled, prec0); break; case 1: - acb_hypgeom_bessel_k_0f1(w0, nu0, z, prec0); + acb_hypgeom_bessel_k_0f1(w0, nu0, z, scaled, prec0); break; default: - acb_hypgeom_bessel_k(w0, nu0, z, prec0); + if (scaled) + acb_hypgeom_bessel_k_scaled(w0, nu0, z, prec0); + else + acb_hypgeom_bessel_k(w0, nu0, z, prec0); } switch (n_randint(state, 3)) { case 0: - acb_hypgeom_bessel_k_asymp(w1, nu0, z, prec1); + acb_hypgeom_bessel_k_asymp(w1, nu0, z, scaled, prec1); break; case 1: - acb_hypgeom_bessel_k_0f1(w1, nu0, z, prec1); + acb_hypgeom_bessel_k_0f1(w1, nu0, z, scaled, prec1); break; default: - acb_hypgeom_bessel_k(w1, nu0, z, prec1); + if (scaled) + acb_hypgeom_bessel_k_scaled(w1, nu0, z, prec1); + else + acb_hypgeom_bessel_k(w1, nu0, z, prec1); } if (!acb_overlaps(w0, w1)) { flint_printf("FAIL: consistency\n\n"); + flint_printf("scaled = %d\n\n", scaled); flint_printf("nu = "); acb_printd(nu0, 30); flint_printf("\n\n"); flint_printf("z = "); acb_printd(z, 30); flint_printf("\n\n"); flint_printf("w0 = "); acb_printd(w0, 30); flint_printf("\n\n"); @@ -89,25 +98,40 @@ int main() switch (n_randint(state, 3)) { case 0: - acb_hypgeom_bessel_k_asymp(w1, nu1, z, prec1); + acb_hypgeom_bessel_k_asymp(w1, nu1, z, scaled, prec1); break; case 1: - acb_hypgeom_bessel_k_0f1(w1, nu1, z, prec1); + acb_hypgeom_bessel_k_0f1(w1, nu1, z, scaled, prec1); break; default: - acb_hypgeom_bessel_k(w1, nu1, z, prec1); + if (scaled) + acb_hypgeom_bessel_k_scaled(w1, nu1, z, prec1); + else + acb_hypgeom_bessel_k(w1, nu1, z, prec1); } switch (n_randint(state, 3)) { case 0: - acb_hypgeom_bessel_k_asymp(w2, nu2, z, prec2); + acb_hypgeom_bessel_k_asymp(w2, nu2, z, scaled, prec2); break; case 1: - acb_hypgeom_bessel_k_0f1(w2, nu2, z, prec2); + acb_hypgeom_bessel_k_0f1(w2, nu2, z, scaled, prec2); break; default: - acb_hypgeom_bessel_k(w2, nu2, z, prec2); + if (scaled) + { + acb_hypgeom_bessel_k(w2, nu2, z, prec2); + acb_exp(t, z, prec2); + acb_mul(w2, w2, t, prec2); + } + else + { + acb_hypgeom_bessel_k_scaled(w2, nu2, z, prec2); + acb_neg(t, z); + acb_exp(t, t, prec2); + acb_mul(w2, w2, t, prec2); + } } acb_mul(t, w1, nu1, prec0); diff --git a/arb_hypgeom.h b/arb_hypgeom.h index 28598059..cadc95bd 100644 --- a/arb_hypgeom.h +++ b/arb_hypgeom.h @@ -74,6 +74,9 @@ void arb_hypgeom_bessel_jy(arb_t res1, arb_t res2, const arb_t nu, const arb_t z void arb_hypgeom_bessel_i(arb_t res, const arb_t nu, const arb_t z, slong prec); void arb_hypgeom_bessel_k(arb_t res, const arb_t nu, const arb_t z, slong prec); +void arb_hypgeom_bessel_i_scaled(arb_t res, const arb_t nu, const arb_t z, slong prec); +void arb_hypgeom_bessel_k_scaled(arb_t res, const arb_t nu, const arb_t z, slong prec); + void arb_hypgeom_airy(arb_t ai, arb_t aip, arb_t bi, arb_t bip, const arb_t z, slong prec); void arb_hypgeom_airy_jet(arb_ptr ai, arb_ptr bi, const arb_t z, slong len, slong prec); void arb_hypgeom_airy_series(arb_poly_t ai, arb_poly_t ai_prime, diff --git a/arb_hypgeom/test/t-wrappers.c b/arb_hypgeom/test/t-wrappers.c index 717a909c..09379987 100644 --- a/arb_hypgeom/test/t-wrappers.c +++ b/arb_hypgeom/test/t-wrappers.c @@ -194,10 +194,18 @@ int main() arb_set_str(v, "[0.75761498638991 +/- 4.63e-15]", prec); TEST(r, v, "bessel_i"); + arb_hypgeom_bessel_i_scaled(r, a, z, prec); + arb_set_str(v, "[0.357871979425934 +/- 8.04e-16]", prec); + TEST(r, v, "bessel_i"); + arb_hypgeom_bessel_k(r, a, z, prec); arb_set_str(v, "[0.68361006034952 +/- 7.89e-15]", prec); TEST(r, v, "bessel_k"); + arb_hypgeom_bessel_k_scaled(r, a, z, prec); + arb_set_str(v, "[1.4472025091165 +/- 4.23e-14]", prec); + TEST(r, v, "bessel_k"); + arb_hypgeom_airy(r, NULL, NULL, NULL, z, prec); arb_set_str(v, "[0.179336305478645 +/- 2.36e-16]", prec); TEST(r, v, "airy ai"); diff --git a/arb_hypgeom/wrappers.c b/arb_hypgeom/wrappers.c index ec3714ca..051312da 100644 --- a/arb_hypgeom/wrappers.c +++ b/arb_hypgeom/wrappers.c @@ -373,6 +373,23 @@ arb_hypgeom_bessel_i(arb_t res, const arb_t nu, const arb_t z, slong prec) acb_clear(u); } +void +arb_hypgeom_bessel_i_scaled(arb_t res, const arb_t nu, const arb_t z, slong prec) +{ + acb_t t, u; + acb_init(t); + acb_init(u); + arb_set(acb_realref(t), nu); + arb_set(acb_realref(u), z); + acb_hypgeom_bessel_i_scaled(t, t, u, prec); + if (acb_is_finite(t) && acb_is_real(t)) + arb_swap(res, acb_realref(t)); + else + arb_indeterminate(res); + acb_clear(t); + acb_clear(u); +} + void arb_hypgeom_bessel_k(arb_t res, const arb_t nu, const arb_t z, slong prec) { @@ -390,6 +407,23 @@ arb_hypgeom_bessel_k(arb_t res, const arb_t nu, const arb_t z, slong prec) acb_clear(u); } +void +arb_hypgeom_bessel_k_scaled(arb_t res, const arb_t nu, const arb_t z, slong prec) +{ + acb_t t, u; + acb_init(t); + acb_init(u); + arb_set(acb_realref(t), nu); + arb_set(acb_realref(u), z); + acb_hypgeom_bessel_k_scaled(t, t, u, prec); + if (acb_is_finite(t) && acb_is_real(t)) + arb_swap(res, acb_realref(t)); + else + arb_indeterminate(res); + acb_clear(t); + acb_clear(u); +} + void arb_hypgeom_expint(arb_t res, const arb_t s, const arb_t z, slong prec) { diff --git a/doc/source/acb_hypgeom.rst b/doc/source/acb_hypgeom.rst index d6aceac3..e82b9a0d 100644 --- a/doc/source/acb_hypgeom.rst +++ b/doc/source/acb_hypgeom.rst @@ -413,12 +413,17 @@ Bessel functions construct the Bessel functions of the third kind (Hankel functions) `H_{\nu}^{(1)}(z), H_{\nu}^{(2)}(z) = J_{\nu}(z) \pm i Y_{\nu}(z)`. -.. function:: void acb_hypgeom_bessel_i_asymp(acb_t res, const acb_t nu, const acb_t z, slong prec) +Modified Bessel functions +------------------------------------------------------------------------------- -.. function:: void acb_hypgeom_bessel_i_0f1(acb_t res, const acb_t nu, const acb_t z, slong prec) +.. function:: void acb_hypgeom_bessel_i_asymp(acb_t res, const acb_t nu, const acb_t z, int scaled, slong prec) + +.. function:: void acb_hypgeom_bessel_i_0f1(acb_t res, const acb_t nu, const acb_t z, int scaled, slong prec) .. function:: void acb_hypgeom_bessel_i(acb_t res, const acb_t nu, const acb_t z, slong prec) +.. function:: void acb_hypgeom_bessel_i_scaled(acb_t res, const acb_t nu, const acb_t z, slong prec) + Computes the modified Bessel function of the first kind `I_{\nu}(z) = z^{\nu} (iz)^{-\nu} J_{\nu}(iz)` respectively using asymptotic series (see :func:`acb_hypgeom_bessel_j_asymp`), @@ -431,7 +436,10 @@ Bessel functions or an automatic algorithm choice. -.. function:: void acb_hypgeom_bessel_k_asymp(acb_t res, const acb_t nu, const acb_t z, slong prec) + The *scaled* version computes the function `e^{-z} I_{\nu}(z)`. The *asymp* + and *0f1* functions implement both variants and allow choosing with a flag. + +.. function:: void acb_hypgeom_bessel_k_asymp(acb_t res, const acb_t nu, const acb_t z, int scaled, slong prec) Computes the modified Bessel function of the second kind via via :func:`acb_hypgeom_u_asymp`. For all `\nu` and all `z \ne 0`, we have @@ -441,7 +449,9 @@ Bessel functions K_{\nu}(z) = \left(\frac{2z}{\pi}\right)^{-1/2} e^{-z} U^{*}(\nu+\tfrac{1}{2}, 2\nu+1, 2z). -.. function:: void acb_hypgeom_bessel_k_0f1_series(acb_poly_t res, const acb_poly_t nu, const acb_poly_t z, slong len, slong prec) + If *scaled* is set, computes the function `e^{z} K_{\nu}(z)`. + +.. function:: void acb_hypgeom_bessel_k_0f1_series(acb_poly_t res, const acb_poly_t nu, const acb_poly_t z, int scaled, slong len, slong prec) Computes the modified Bessel function of the second kind `K_{\nu}(z)` as a power series truncated to length *len*, @@ -462,7 +472,9 @@ Bessel functions As currently implemented, the output is indeterminate if `\nu[0]` is nonexact and contains an integer. -.. function:: void acb_hypgeom_bessel_k_0f1(acb_t res, const acb_t nu, const acb_t z, slong prec) + If *scaled* is set, computes the function `e^{z} K_{\nu}(z)`. + +.. function:: void acb_hypgeom_bessel_k_0f1(acb_t res, const acb_t nu, const acb_t z, int scaled, slong prec) Computes the modified Bessel function of the second kind from @@ -483,11 +495,17 @@ Bessel functions As currently implemented, the output is indeterminate if `\nu` is nonexact and contains an integer. + If *scaled* is set, computes the function `e^{z} K_{\nu}(z)`. + .. function:: void acb_hypgeom_bessel_k(acb_t res, const acb_t nu, const acb_t z, slong prec) Computes the modified Bessel function of the second kind `K_{\nu}(z)` using an automatic algorithm choice. +.. function:: void acb_hypgeom_bessel_k_scaled(acb_t res, const acb_t nu, const acb_t z, slong prec) + + Computes the function `e^{z} K_{\nu}(z)`. + Airy functions ------------------------------------------------------------------------------- diff --git a/doc/source/arb_hypgeom.rst b/doc/source/arb_hypgeom.rst index e4b4283f..8d7da401 100644 --- a/doc/source/arb_hypgeom.rst +++ b/doc/source/arb_hypgeom.rst @@ -278,10 +278,18 @@ Bessel functions Computes the modified Bessel function of the first kind `I_{\nu}(z) = z^{\nu} (iz)^{-\nu} J_{\nu}(iz)`. +.. function:: void arb_hypgeom_bessel_i_scaled(arb_t res, const arb_t nu, const arb_t z, slong prec) + + Computes the function `e^{-z} I_{\nu}(z)`. + .. function:: void arb_hypgeom_bessel_k(arb_t res, const arb_t nu, const arb_t z, slong prec) Computes the modified Bessel function of the second kind `K_{\nu}(z)`. +.. function:: void arb_hypgeom_bessel_k_scaled(arb_t res, const arb_t nu, const arb_t z, slong prec) + + Computes the function `e^{z} K_{\nu}(z)`. + Airy functions -------------------------------------------------------------------------------