ENH: fix j1 in the far left tail

[dschmitz89]

Reference issue

Towards scipy/scipy#22706

What does this implement/fix?

Fixes the $j1$ function in the far left tail using the suggestion $x/2$ for $x<sqrt(eps)$.

Additional information

I attach a crude testing script that establishes $sqrt(epsilon)$ as a suitable cutoff point.

Details

import numpy as np
import matplotlib.pyplot as plt
from mpmath import mp
from scipy.special import j1

mp.dps = 1000

@np.vectorize
def mpmath_j1(x):
    x = mp.mpf(x)
    return float(mp.besselj(1, x))

min_subnormal = 2**(-1074)
x_values = np.logspace(np.log10(min_subnormal), 1, num=5000)
mpmath_results = mpmath_j1(x_values)
scipy_results = j1(x_values)
x_half = 0.5 * x_values
relative_errors_j1 = np.abs(mpmath_results - scipy_results) / np.abs(mpmath_results)
relative_errors_x_half = np.abs(mpmath_results - x_half) / np.abs(mpmath_results)
j1_better = relative_errors_j1 < relative_errors_x_half
print(f"First x value where j1 is better: {x_values[j1_better][0]}")
plt.figure(figsize=(10, 6))
#plt.loglog(x_values, relative_errors_j1, label='j1', color='blue')
plt.loglog(x_values, relative_errors_x_half, label='0.5 * x', color='green', linestyle='dashed')
#plt.title('Comparison of mpmath.j1 and scipy.special.j1')
plt.xlabel('x')
plt.ylabel('relative error')
plt.grid()
plt.legend()
plt.show()

AI Generation Disclosure

AI was used to write the test cases. The actual test case values were manually created.

[fbourgey]

I did a quick check and this looks good to me:

• The PR is a bit of a misnomer. I would avoid far left tail and use something like ENH: use leading-order approximation for j1 at small |x|.
• I guess the same patch could be used for j0, y1, etc?