diff --git a/doc/source/acb_hypgeom.rst b/doc/source/acb_hypgeom.rst
index e2c5c25375..13299b7114 100644
--- a/doc/source/acb_hypgeom.rst
+++ b/doc/source/acb_hypgeom.rst
@@ -445,6 +445,11 @@ Error functions and Fresnel integrals
Bessel functions
-------------------------------------------------------------------------------
+.. function:: void acb_hypgeom_bessel_j_deriv_bound(mag_t res, const acb_t nu, const acb_t z, ulong d)
+
+ Sets *res* to a bound, possibly crude, for `|J^{(d)}_{\nu}(z)|`.
+ Currently only specialized for small integer `\nu` and small `d`.
+
.. function:: void acb_hypgeom_bessel_j_asymp(acb_t res, const acb_t nu, const acb_t z, slong prec)
Computes the Bessel function of the first kind
diff --git a/examples/integrals.c b/examples/integrals.c
index d93d3ae580..a4da4f5eaf 100644
--- a/examples/integrals.c
+++ b/examples/integrals.c
@@ -743,11 +743,33 @@ f_si(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
return 0;
}
+/* f(z) = J_0(z) J_1(z) J_2(z) */
+int
+f_triple_bessel_j(acb_ptr res, const acb_t z, void * param, slong order, slong prec)
+{
+ acb_t t;
+
+ if (order > 1)
+ flint_abort(); /* Would be needed for Taylor method. */
+
+ acb_init(t);
+ acb_hypgeom_bessel_j(res, t, z, prec);
+ acb_set_ui(t, 1);
+ acb_hypgeom_bessel_j(t, t, z, prec);
+ acb_mul(res, res, t, prec);
+ acb_set_ui(t, 2);
+ acb_hypgeom_bessel_j(t, t, z, prec);
+ acb_mul(res, res, t, prec);
+ acb_clear(t);
+
+ return 0;
+}
+
/* ------------------------------------------------------------------------- */
/* Main test program */
/* ------------------------------------------------------------------------- */
-#define NUM_INTEGRALS 38
+#define NUM_INTEGRALS 39
const char * descr[NUM_INTEGRALS] =
{
@@ -789,6 +811,7 @@ const char * descr[NUM_INTEGRALS] =
"int_{-1-i}^{-1+i} 1/sqrt(x) dx",
"int_0^{inf} 1/gamma(x) dx (using domain truncation)",
"int_0^{1000} Si(x) dx",
+ "int_0^{1000} J_0(x) J_1(x) J_2(x) dx",
};
int main(int argc, char *argv[])
@@ -1291,6 +1314,11 @@ int main(int argc, char *argv[])
acb_calc_integrate(s, f_si, NULL, a, b, goal, tol, options, prec);
break;
+ case 38:
+ acb_zero(a);
+ acb_set_ui(b, 1000);
+ acb_calc_integrate(s, f_triple_bessel_j, NULL, a, b, goal, tol, options, prec);
+ break;
default:
abort();
diff --git a/src/acb_hypgeom.h b/src/acb_hypgeom.h
index d87595842b..9d7ca66eeb 100644
--- a/src/acb_hypgeom.h
+++ b/src/acb_hypgeom.h
@@ -130,6 +130,7 @@ void acb_hypgeom_m_1f1(acb_t res, const acb_t a, const acb_t b, const acb_t z, i
void acb_hypgeom_m(acb_t res, const acb_t a, const acb_t b, const acb_t z, int regularized, slong prec);
void acb_hypgeom_1f1(acb_t res, const acb_t a, const acb_t b, const acb_t z, int regularized, slong prec);
+void acb_hypgeom_bessel_j_deriv_bound(mag_t res, const acb_t nu, const acb_t z, ulong d);
void acb_hypgeom_bessel_j_0f1(acb_t res, const acb_t nu, const acb_t z, slong prec);
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);
diff --git a/src/acb_hypgeom/bessel_j.c b/src/acb_hypgeom/bessel_j.c
index e19bdba769..2c56de2499 100644
--- a/src/acb_hypgeom/bessel_j.c
+++ b/src/acb_hypgeom/bessel_j.c
@@ -1,5 +1,5 @@
/*
- Copyright (C) 2014-2015 Fredrik Johansson
+ Copyright (C) 2014-2015, 2025 Fredrik Johansson
This file is part of FLINT.
@@ -9,9 +9,253 @@
(at your option) any later version. See .
*/
+#include
#include "acb.h"
#include "acb_hypgeom.h"
+/*
+Bound |J^{(d)}_nu(z)|, or |J_nu(z)| itself with d = 0.
+*/
+void
+acb_hypgeom_bessel_j_deriv_bound(mag_t res, const acb_t nu, const acb_t z, ulong d)
+{
+ if (!acb_is_int(nu))
+ {
+ /* Not yet implemented. */
+ /* Things we could use:
+
+ - For complex z, complex n, integer d >= 0:
+ |J^{(d)}_n(z)| <= max(|J_{n-d}(z)|, ..., |J_{n+d}(z)|).
+ - Asymptotic expansion (general case)
+ Bounds are messy and not useful when |nu| is large.
+ */
+ mag_inf(res);
+ return;
+ }
+
+ mag_t bound, t, u, x, y, zlow;
+
+ mag_init(bound);
+ mag_init(t);
+ mag_init(u);
+ mag_init(x);
+ mag_init(y);
+ mag_init(zlow);
+
+ arb_get_mag(x, acb_realref(z));
+ arb_get_mag(y, acb_imagref(z));
+ acb_get_mag_lower(zlow, z);
+
+ /* |J^{(d)}_n(x)| <= exp(|im(z)|). */
+ mag_exp(bound, y);
+
+ /* |n| is something reasonable, do something less generic */
+ if (arf_cmpabs_2exp_si(arb_midref(acb_realref(nu)), 30) < 0)
+ {
+ slong n;
+
+ n = arf_get_si(arb_midref(acb_realref(nu)), ARF_RND_DOWN);
+ n = FLINT_ABS(n);
+
+ /*
+ Small |z| bound.
+ [DLMF 10.14.4] For real n >= -0.5, complex z:
+
+ |J_n(z)| <= (|z/2|^n / n!) exp(|im(z)|)
+
+ Analogously for d = 1:
+
+ |J'_0(z)| <= (|z|/2) exp(|im(z)|)
+ |J'_1(z)| <= (1/2) exp(|im(z)|)
+ |J'_n(z)| <= (|z|^(n-1) / (2^n * (n-1)!)) exp(|im(z)|)
+
+ TODO: generalize to higher derivatives.
+ */
+ if (d == 0)
+ {
+ mag_hypot(t, x, y);
+ mag_mul_2exp_si(t, t, -1);
+ mag_pow_ui(t, t, n);
+ mag_rfac_ui(u, n);
+ mag_mul(t, t, u);
+ if (mag_cmp_2exp_si(t, 0) < 0)
+ mag_mul(bound, bound, t);
+ }
+ else if (d == 1)
+ {
+ if (n == 1)
+ {
+ mag_mul_2exp_si(bound, bound, -1);
+ }
+ else
+ {
+ mag_hypot(t, x, y);
+
+ if (n == 0)
+ mag_mul_2exp_si(t, t, -1);
+ else
+ {
+ mag_pow_ui(t, t, n - 1);
+ mag_mul_2exp_si(t, t, -n);
+ mag_rfac_ui(u, n - 1);
+ mag_mul(t, t, u);
+ }
+
+ if (mag_cmp_2exp_si(t, 0) < 0)
+ mag_mul(bound, bound, t);
+ }
+ }
+
+ /*
+ For integer n,
+
+ |J^{(d)}_n(z)| <= C_{d,n} * cosh(im(z)) / sqrt(2/(pi |z|))
+
+ for some C_{d,n} >= 1.
+
+ Proving tight explicit bounds for C_{d,n} is an open problem:
+ https://mathoverflow.net/questions/485404
+
+ However, we can brute force some reasonable bounds for small n, d
+ which are most commonly needed. TODO: extend this table.
+ */
+ if (n <= 16 && d <= 16)
+ {
+ /* sqrt(2/pi) + 0.0001 */
+ double Csq2pi = 0.798;
+
+ static const float Cbound[] = {
+ 1.077, /* n = 0 */
+ 1.035, /* n = 1 */
+ 1.089, /* n = 2 */
+ 1.131, /* n = 3 */
+ 1.167, /* n = 4 */
+ 1.197, /* n = 5 */
+ 1.223, /* n = 6 */
+ 1.247, /* n = 7 */
+ 1.269, /* n = 8 */
+ 1.288, /* n = 9 */
+ 1.306, /* n = 10 */
+ 1.323, /* n = 11 */
+ 1.339, /* n = 12 */
+ 1.354, /* n = 13 */
+ 1.368, /* n = 14 */
+ 1.382, /* n = 15 */
+ 1.394, /* n = 16 */
+ };
+
+ if (d == 0)
+ Csq2pi *= Cbound[n];
+ else if (n == 0 && d == 1)
+ Csq2pi *= Cbound[1];
+
+ mag_rsqrt(t, zlow);
+ mag_set_d(u, Csq2pi);
+ mag_mul(t, t, u);
+ mag_cosh(u, y);
+ mag_mul(t, t, u);
+ mag_min(bound, bound, t);
+ }
+ else if (mag_is_zero(y))
+ {
+ /* Landau: |J^{(d)}_n(x)| <= 0.786 |x|^(-1/3) */
+ mag_inv(t, zlow);
+ mag_root(t, t, 3);
+ mag_mul_ui(t, t, 13186892);
+ mag_mul_2exp_si(t, t, -24);
+ mag_min(bound, bound, t);
+ }
+
+ /* Landau:
+ |J_n(x)| <= 0.675 |n|^(-1/3)
+ |J^{(d)}_n(x)| <= 0.675 max(0, |n|-d)^(-1/3)
+ */
+ if (mag_is_zero(y) && d < (ulong) n)
+ {
+ mag_set_ui_lower(t, n - d);
+ mag_inv(t, t);
+ mag_root(t, t, 3);
+ mag_mul_ui(t, t, 11324621);
+ mag_mul_2exp_si(t, t, -24);
+ mag_min(bound, bound, t);
+ }
+ }
+
+ mag_set(res, bound);
+
+ mag_clear(bound);
+ mag_clear(t);
+ mag_clear(u);
+ mag_clear(x);
+ mag_clear(y);
+ mag_clear(zlow);
+}
+
+static int
+_acb_hypgeom_bessel_j_is_real(const acb_t nu, const acb_t z)
+{
+ if (acb_is_int(nu))
+ {
+ if (arb_is_zero(acb_imagref(z)))
+ return 1;
+ if (arb_is_zero(acb_realref(z)))
+ return arf_is_int_2exp_si(arb_midref(acb_realref(nu)), 1);
+ }
+
+ return 0;
+}
+
+static int
+_acb_hypgeom_bessel_j_is_imag(const acb_t nu, const acb_t z)
+{
+ if (acb_is_int(nu) && arb_is_zero(acb_realref(z)))
+ return !arf_is_int_2exp_si(arb_midref(acb_realref(nu)), 1);
+ else
+ return 0;
+}
+
+/* Bound propagated error for |J_n(z)| when evaluating at mid(z). */
+void
+_acb_hypgeom_bessel_j_prop_error(mag_t re, mag_t im, const acb_t nu, const acb_t z)
+{
+ mag_t err, rad;
+
+ mag_init(err);
+ mag_init(rad);
+
+ acb_hypgeom_bessel_j_deriv_bound(err, nu, z, 1);
+
+ if (!mag_is_finite(err))
+ {
+ mag_inf(re);
+ mag_inf(im);
+ }
+ else
+ {
+ mag_hypot(rad, arb_radref(acb_realref(z)), arb_radref(acb_imagref(z)));
+ mag_mul(err, err, rad);
+
+ if (_acb_hypgeom_bessel_j_is_real(nu, z))
+ {
+ mag_set(re, err);
+ mag_zero(im);
+ }
+ else if (_acb_hypgeom_bessel_j_is_imag(nu, z))
+ {
+ mag_zero(re);
+ mag_set(im, err);
+ }
+ else
+ {
+ mag_set(re, err);
+ mag_set(im, err);
+ }
+ }
+
+ mag_clear(err);
+ mag_clear(rad);
+}
+
/* assumes no aliasing */
/* (+/- iz)^(-1/2-v) * z^v * exp(+/- iz) */
void
@@ -175,8 +419,8 @@ acb_hypgeom_bessel_j_asymp(acb_t res, const acb_t nu, const acb_t z, slong prec)
acb_clear(u);
}
-void
-acb_hypgeom_bessel_j_0f1(acb_t res, const acb_t nu, const acb_t z, slong prec)
+static void
+_acb_hypgeom_bessel_j_0f1(acb_t res, const acb_t nu, const acb_t z, slong prec, slong prec2)
{
acb_struct b[2];
acb_t w, c, t;
@@ -186,7 +430,7 @@ acb_hypgeom_bessel_j_0f1(acb_t res, const acb_t nu, const acb_t z, slong prec)
acb_init(t);
acb_neg(t, nu);
- acb_hypgeom_bessel_j_0f1(res, t, z, prec);
+ _acb_hypgeom_bessel_j_0f1(res, t, z, prec, prec2);
acb_mul_2exp_si(t, t, -1);
if (!acb_is_int(t))
@@ -202,7 +446,7 @@ acb_hypgeom_bessel_j_0f1(acb_t res, const acb_t nu, const acb_t z, slong prec)
acb_init(c);
acb_init(t);
- acb_add_ui(b + 0, nu, 1, prec);
+ acb_add_ui(b + 0, nu, 1, prec2);
acb_one(b + 1);
/* (z/2)^nu / gamma(nu+1) */
@@ -212,11 +456,11 @@ acb_hypgeom_bessel_j_0f1(acb_t res, const acb_t nu, const acb_t z, slong prec)
acb_mul(c, t, c, prec);
/* -z^2/4 */
- acb_mul(w, z, z, prec);
+ acb_mul(w, z, z, prec2);
acb_mul_2exp_si(w, w, -2);
acb_neg(w, w);
- acb_hypgeom_pfq_direct(t, NULL, 0, b, 2, w, -1, prec);
+ acb_hypgeom_pfq_direct(t, NULL, 0, b, 2, w, -1, prec2);
acb_mul(res, t, c, prec);
@@ -227,28 +471,180 @@ acb_hypgeom_bessel_j_0f1(acb_t res, const acb_t nu, const acb_t z, slong prec)
acb_clear(t);
}
-/*
+void
+acb_hypgeom_bessel_j_0f1(acb_t res, const acb_t nu, const acb_t z, slong prec)
+{
+ _acb_hypgeom_bessel_j_0f1(res, nu, z, prec, prec);
+}
+
+void
+acb_hypgeom_bessel_j_transition(acb_t res, const acb_t nu, const acb_t z, slong prec)
+{
+ slong prec2;
+ double a, b, zz;
+ slong cancellation;
+
+ /* All terms are positive, so no cancellation */
+ if (arb_is_zero(acb_realref(z)) && (acb_is_int(nu) || (acb_is_real(nu) && arb_is_positive(acb_realref(nu)))))
+ {
+ acb_hypgeom_bessel_j_0f1(res, nu, z, prec);
+ return;
+ }
-The asymptotic series can be used roughly when
+ a = arf_get_d(arb_midref(acb_realref(z)), ARF_RND_DOWN);
+ b = arf_get_d(arb_midref(acb_imagref(z)), ARF_RND_DOWN);
+ zz = sqrt(a * a + b * b);
-[(1+log(2))/log(2) = 2.44269504088896] * z > p
+ /* If nu > |z|^2/4, there is no significant cancellation. */
+ if (acb_is_real(nu) && arf_cmpabs_2exp_si(arb_midref(acb_realref(nu)), 20) < 64)
+ {
+ double nn;
+ nn = arf_get_d(arb_midref(acb_realref(nu)), ARF_RND_DOWN);
+
+ if (nn > zz * zz * 0.25)
+ {
+ acb_hypgeom_bessel_j_0f1(res, nu, z, prec);
+ return;
+ }
+ }
-We are a bit more conservative and use the factor 2.
-*/
+ /* Estimate cancellation as (|x|-im(|x|)) * (1/log(2)) bits.
+ TODO: estimate more accurately for large nu */
+ cancellation = (zz - fabs(b)) * 1.44269504088896;
+ cancellation = FLINT_MIN(cancellation, 4 * prec);
+ cancellation = FLINT_MAX(cancellation, 0);
+
+ prec2 = prec + 5 + cancellation;
+
+ if (acb_is_exact(nu) && acb_is_exact(z))
+ {
+ _acb_hypgeom_bessel_j_0f1(res, nu, z, prec, prec2);
+ }
+ else
+ {
+ mag_t aerr, berr;
+ acb_t t;
+
+ mag_init(aerr);
+ mag_init(berr);
+
+ _acb_hypgeom_bessel_j_prop_error(aerr, berr, nu, z);
+
+ /* Propagated error implemented for this nu. */
+ if (mag_is_finite(aerr) && mag_is_finite(berr))
+ {
+ acb_init(t);
+ acb_get_mid(t, z);
+
+ _acb_hypgeom_bessel_j_0f1(res, nu, t, prec, prec2);
+
+ arb_add_error_mag(acb_realref(res), aerr);
+ arb_add_error_mag(acb_imagref(res), berr);
+ acb_clear(t);
+ }
+ else /* Propagated error not implemented */
+ {
+ slong acc = acb_rel_accuracy_bits(z);
+
+ prec2 = FLINT_MIN(prec2, acc);
+ prec2 = FLINT_MAX(prec2, prec);
+
+ _acb_hypgeom_bessel_j_0f1(res, nu, z, prec, prec2);
+ }
+
+ mag_clear(aerr);
+ mag_clear(berr);
+ }
+}
void
acb_hypgeom_bessel_j(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_j_0f1(res, nu, z, prec);
+ if (!acb_is_finite(nu) || !acb_is_finite(z))
+ {
+ /* Some infinities */
+ if (acb_is_real(nu) && acb_is_finite(nu) &&
+ acb_is_real(z) && (mag_is_finite(arb_radref(acb_realref(z))) && arf_is_inf(arb_midref(acb_realref(z)))))
+ acb_zero(res);
+ else
+ acb_indeterminate(res);
+ }
else
- acb_hypgeom_bessel_j_asymp(res, nu, z, prec);
+ {
+ mag_t zmag;
+
+ mag_init(zmag);
+ acb_get_mag(zmag, z);
- mag_clear(zmag);
+ if (mag_cmp_2exp_si(zmag, 3) < 0)
+ {
+ acb_hypgeom_bessel_j_0f1(res, nu, z, prec);
+ }
+ else
+ {
+ acb_t t;
+ acb_init(t);
+
+ /* Assuming small nu, the asymptotic series can be used roughly when
+ [(1+log(2))/log(2) = 2.44269504088896] * z > p
+ We are a bit more conservative and use the factor 2. */
+
+ if (mag_cmp_2exp_si(zmag, 64) >= 0 || 2 * mag_get_d(zmag) >= prec)
+ acb_hypgeom_bessel_j_asymp(t, nu, z, prec);
+ else
+ acb_hypgeom_bessel_j_transition(t, nu, z, prec);
+
+#if 0
+ flint_printf("T: "); acb_printd(t, 10); flint_printf("\n");
+#endif
+
+ /* If the enclosure is terrible, try bounding |J_nu(z)|
+ directly. TODO: cleaner way to do this (detect
+ wide intervals a priori?) */
+ if (acb_rel_accuracy_bits(t) < 1)
+ {
+ mag_t M;
+ int real, imag;
+
+ mag_init(M);
+
+ acb_hypgeom_bessel_j_deriv_bound(M, nu, z, 0);
+
+#if 0
+ flint_printf("M: "); mag_printd(M, 10); flint_printf("\n");
+#endif
+
+ if (mag_is_finite(M))
+ {
+ real = _acb_hypgeom_bessel_j_is_real(nu, z);
+ imag = _acb_hypgeom_bessel_j_is_imag(nu, z);
+
+ acb_zero(res);
+ if (real)
+ arb_add_error_mag(acb_realref(res), M);
+ else if (imag)
+ arb_add_error_mag(acb_imagref(res), M);
+ else
+ acb_add_error_mag(res, M);
+
+ if (acb_is_finite(t))
+ {
+ arb_intersection(acb_realref(t), acb_realref(t), acb_realref(res), prec);
+ arb_intersection(acb_imagref(t), acb_imagref(t), acb_imagref(res), prec);
+ }
+ else
+ {
+ acb_set(t, res);
+ }
+ }
+
+ mag_clear(M);
+ }
+
+ acb_swap(res, t);
+ acb_clear(t);
+ }
+
+ mag_clear(zmag);
+ }
}
diff --git a/src/acb_hypgeom/test/main.c b/src/acb_hypgeom/test/main.c
index 176790f8f2..20bbd48234 100644
--- a/src/acb_hypgeom/test/main.c
+++ b/src/acb_hypgeom/test/main.c
@@ -19,6 +19,7 @@
#include "t-airy_series.c"
#include "t-bessel_i.c"
#include "t-bessel_j.c"
+#include "t-bessel_j_deriv_bound.c"
#include "t-bessel_k.c"
#include "t-bessel_y.c"
#include "t-beta_lower.c"
@@ -87,6 +88,7 @@ test_struct tests[] =
TEST_FUNCTION(acb_hypgeom_airy_series),
TEST_FUNCTION(acb_hypgeom_bessel_i),
TEST_FUNCTION(acb_hypgeom_bessel_j),
+ TEST_FUNCTION(acb_hypgeom_bessel_j_deriv_bound),
TEST_FUNCTION(acb_hypgeom_bessel_k),
TEST_FUNCTION(acb_hypgeom_bessel_y),
TEST_FUNCTION(acb_hypgeom_beta_lower),
diff --git a/src/acb_hypgeom/test/t-bessel_j_deriv_bound.c b/src/acb_hypgeom/test/t-bessel_j_deriv_bound.c
new file mode 100644
index 0000000000..551764d5c1
--- /dev/null
+++ b/src/acb_hypgeom/test/t-bessel_j_deriv_bound.c
@@ -0,0 +1,147 @@
+/*
+ Copyright (C) 2025 Fredrik Johansson
+
+ This file is part of FLINT.
+
+ FLINT is free software: you can redistribute it and/or modify it under
+ the terms of the GNU Lesser General Public License (LGPL) as published
+ by the Free Software Foundation; either version 3 of the License, or
+ (at your option) any later version. See .
+*/
+
+#include "test_helpers.h"
+#include "acb.h"
+#include "acb_hypgeom.h"
+
+static void
+_acb_randtest_inner(acb_t z2, flint_rand_t state, const acb_t z1)
+{
+ 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)));
+ }
+
+ if (!acb_contains(z1, z2))
+ flint_abort();
+}
+
+TEST_FUNCTION_START(acb_hypgeom_bessel_j_deriv_bound, state)
+{
+ slong iter;
+
+ for (iter = 0; iter < 10000 * 0.1 * flint_test_multiplier(); iter++)
+ {
+ acb_t f, nu1, nu2, z1, z2;
+ slong prec;
+ mag_t B, fB;
+ ulong d;
+
+ acb_init(z1);
+ acb_init(z2);
+ acb_init(nu1);
+ acb_init(nu2);
+ acb_init(f);
+ mag_init(B);
+ mag_init(fB);
+
+ 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_randtest_inner(z2, state, z1);
+
+ /* Currently only integer nu is implemented, so don't
+ waste time testing anything else. */
+ if (n_randint(state, 2))
+ acb_set_si(nu1, n_randint(state, 10));
+ else
+ acb_set_si(nu1, n_randint(state, 100));
+
+ if (n_randint(state, 2))
+ acb_neg(nu1, nu1);
+ acb_set(nu2, nu1);
+
+ d = n_randint(state, 2);
+
+ acb_hypgeom_bessel_j_deriv_bound(B, nu1, z1, d);
+
+ prec = MAG_BITS + 10;
+
+ do {
+ if (d == 0)
+ {
+ acb_hypgeom_bessel_j(f, nu2, z2, prec);
+ }
+ else
+ {
+ /* J'_nu = (J_{nu-1} - J_{nu+1})/2 */
+ acb_t t;
+ acb_init(t);
+ acb_sub_ui(t, nu2, 1, prec);
+ acb_hypgeom_bessel_j(f, t, z2, prec);
+ acb_add_ui(t, nu2, 1, prec);
+ acb_hypgeom_bessel_j(t, t, z2, prec);
+ acb_sub(f, f, t, prec);
+ acb_mul_2exp_si(f, f, -1);
+ acb_clear(t);
+ }
+
+ if (acb_rel_accuracy_bits(f) >= MAG_BITS)
+ break;
+
+ prec *= 2;
+ } while (1);
+
+ acb_get_mag_lower(fB, f);
+
+ if (mag_cmp(fB, B) > 0)
+ {
+ printf("FAIL\n");
+
+ flint_printf("d = %wu\n\n", d);
+ flint_printf("z1 = "); acb_printd(z1, 20); flint_printf("\n");
+ flint_printf("z2 = "); acb_printd(z2, 20); flint_printf("\n\n");
+ flint_printf("nu1 = "); acb_printd(nu1, 20); flint_printf("\n");
+ flint_printf("nu2 = "); acb_printd(nu2, 20); flint_printf("\n\n");
+
+ flint_printf("B = "); mag_printd(B, 10); printf("\n");
+ flint_printf("f = "); acb_printd(f, 20); flint_printf("\n");
+ flint_printf("fB = "); mag_printd(fB, 10); printf("\n\n");
+
+ flint_abort();
+ }
+
+ acb_clear(z1);
+ acb_clear(z2);
+ acb_clear(nu1);
+ acb_clear(nu2);
+ acb_clear(f);
+ mag_clear(B);
+ mag_clear(fB);
+ }
+
+ TEST_FUNCTION_END(state);
+}
diff --git a/src/python/flint_ctypes.py b/src/python/flint_ctypes.py
index 2cc0ac1698..5fc4abc587 100644
--- a/src/python/flint_ctypes.py
+++ b/src/python/flint_ctypes.py
@@ -1969,6 +1969,10 @@ def bessel_j(ctx, x, y):
"""
>>> RR.bessel_j(2, 3)
[0.486091260585891 +/- 4.75e-16]
+ >>> sum(CC.bessel_j(0, k) for k in range(100))
+ [1.419207859380 +/- 2.35e-13]
+ >>> w = CC.exp_pi_i(QQ(1)/100); sum(CC.bessel_j(0, w*k) for k in range(101))
+ ([0.78030446659 +/- 2.44e-12] + [0.62756344686 +/- 5.49e-12]*I)
"""
return ctx._binary_op(x, y, libgr.gr_bessel_j, "bessel_j($n, $x)")