Merge pull request flintlib#2137 from fredrik-johansson/bessel

Improve bessel_j in transition region
albinahlback · Jan 7, 2025 · 151f514 · 151f514
2 parents 86a8c10 + 986bb94
commit 151f514
Show file tree

Hide file tree

Showing 7 changed files with 606 additions and 23 deletions.
diff --git 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
@@ -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
@@ -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);