Product integrals by Mellin residue summation

\(\int_0^\infty e^{-x} J_0(w x)\,dx = 1/\sqrt{1+w^2}\). For large w the integrand oscillates many times per unit length, while the residue series converges like \((1/w)^n\). The example also shows a logarithmic case (coinciding poles), a three-factor product with no closed form, and the engine refusing a case at its convergence boundary instead of guessing.

w =       10: rel. error      0.0     2.7 ms  (41 terms, expand exp(-1 x))
w =   1000.0: rel. error      0.0     1.5 ms  (21 terms, expand exp(-1 x))
w = 100000.0: rel. error  1.5e-31     1.2 ms  (17 terms, expand exp(-1 x))
Int e^(-2x)/(x+3) dx = 0.145267629233886893806035404752  exact: 0.145267629233886893806035404752
Int e^(-2x) e^(-x^2/2) J1(100x) dx = 0.00980001000675268801870689780916
k = b: neither ordering of exp(-1 x) * J_0(1 x) is a convergent residue expansion (the scales are ...

import time

import mpmath as mp

from ramanujan_qft import (SeriesConvergenceError, bessel_j, exponential, gaussian_bessel,
                           inverse_power, product_integral)

mp.mp.dps = 30
for w in (10, 1e3, 1e5):
    t = time.perf_counter()
    v, info = product_integral(exponential(1), bessel_j(0, w), return_info=True)
    dt = time.perf_counter() - t
    err = abs(v - 1 / mp.sqrt(1 + mp.mpf(w)**2)) / v
    print(f"w = {w:>8}: rel. error {mp.nstr(err, 2):>8}  {dt * 1e3:6.1f} ms  ({info['terms']} terms, {info['method']})")

# Logarithmic case (coinciding poles), exact: e^6 E1(6)
print("Int e^(-2x)/(x+3) dx =", product_integral(exponential(2), inverse_power(3, 1)), " exact:", mp.exp(6) * mp.e1(6))

# Three-factor product e^(-2x) e^(-x^2/2) J1(100 x)
print("Int e^(-2x) e^(-x^2/2) J1(100x) dx =", product_integral(exponential(2), gaussian_bessel(0.5, 100, 1)))

# Where the expansion parameter is ~1 the engine refuses instead of guessing:
try:
    product_integral(exponential(1), bessel_j(0, 1), max_terms=300)
except SeriesConvergenceError as e:
    print("k = b:", str(e)[:90], "...")