Static J1_ZEROS 

pub(crate) static J1_ZEROS: [(u64, u64); 48]

J1 zeros and extremums on [-76;76]

Generated in SageMath:

from mpmath import mp, mpf, findroot, j1
from sage.all import *
import struct

DR = RealField(52)
DD = RealField(190)

def double_to_hex(f):
    packed = struct.pack('>d', float(f))
    return '0x' + packed.hex()

def split_double_double(x):
    x_hi = DR(x)  # convert to f64
    x_lo = x - DD(x_hi)
    return (x_lo,x_hi)

def print_double_double(mark, x):
    splat = split_double_double(x)
    print(f"{mark}({double_to_hex(splat[0])}, {double_to_hex(splat[1])}),")

zeros = []

# Step size to detect sign changes
mp.prec = 150

step = mpf("0.001")
epsilon = mpf("1e-35")
x = mpf("0.0")

previous_zero = R120(0)
j1_zeros = []

while x < mpf("76.0"):
    f1 = besselj(1, x)
    f2 = besselj(1, x + step)
    if f1 * f2 < 0:
        zero = findroot(lambda t: j1(t), (x, x + step), solver='bisect', tol=mp.mpf("1e-41"))
        previous_zero = zero
        j1_zeros.append(zero)
    if previous_zero is not None and abs(x - mpf(f'{round(x)}')) < epsilon:
        zeros.append(previous_zero)
    x += step

j1_extrema = []

x = mpf("0.0")
while x < mpf("76.0"):
    d1 = mp.diff(lambda t: j1(t), x)
    d2 = mp.diff(lambda t: j1(t), x + step)
    if d1 * d2 < 0:
        extremum = findroot(lambda t: mp.diff(lambda u: j1(u), t), (x, x + step), solver='bisect', tol=mp.mpf("1e-41"))
        j1_extrema.append(extremum)
    x += step

j1_zeros.extend(j1_extrema)
j1_zeros = sorted(j1_zeros)

print("static J1_ZEROS: [(u64, u64); 46] = [")
for z in j1_zeros:
    k = split_double_double(DD(z))
    hi = double_to_hex(k[1])
    lo = double_to_hex(k[0])
    print(f"({lo}, {hi}),")

print("];")

See notes/bessel_j1_taylor.ipynb