Static J1_COEFFS 

pub(crate) static J1_COEFFS: [[(u64, u64); 24]; 47]

Following search for J1 (see J1_ZEROS) zeros and extremums: at each zero and extremum we’re doing Taylor series expansion one that should be enough to cover whole interval between zero or peak which is PI/4

Generated in SageMath and Sollya:

def compute_intervals(zeros):
    intervals = []
    for i in range(0, len(zeros)):
        if i == 0:
            a = (zeros[i]) / 2 - 0.05 - zeros[i]
            b = (zeros[i] + zeros[i + 1]) / 2 + 0.05 - zeros[i]
            intervals.append((RealField(18)(a), RealField(18)(b), RealField(110)(zeros[i])))
        elif i + 1 > len(zeros) - 1:
            a = (zeros[i - 1] + zeros[i]) / 2 - 0.05 - zeros[i]
            b = (zeros[i]) + 0.83 + 0.05 - zeros[i]
            intervals.append((RealField(18)(a), RealField(18)(b), RealField(110)(zeros[i])))
        else:
            a = (zeros[i - 1] + zeros[i]) / 2 - zeros[i] - 0.05
            b = (zeros[i] + zeros[i + 1]) / 2 + 0.05  - zeros[i]
            intervals.append((RealField(18)(a), RealField(18)(b), RealField(110)(zeros[i])))
    return intervals

intervals = compute_intervals(j1_zeros)
# print(intervals)

def build_sollya_script(a, b, zero, deg):
    return f"""
prec = 500;
bessel_j1 = library("./notes/bessel_sollya/cmake-build-release/libbessel_sollya.dylib");
f = bessel_j1(x + {zero});
d = [{a}, {b}];
pf = remez(f, {deg}, d);
for i from 0 to degree(pf) do {{
    write(coeff(pf, i)) >> "coefficients.txt";
    write("\\n") >> "coefficients.txt";
}};
"""

def load_coefficients(filename):
    with open(filename, "r") as f:
        return [RealField(500)(line.strip()) for line in f if line.strip()]

def call_sollya_on_interval(a, b, zero, degree=12):
    sollya_script = build_sollya_script(a, b, zero, degree)
    with open("tmp_interval.sollya", "w") as f:
        f.write(sollya_script)
    import subprocess
    if os.path.exists("coefficients.txt"):
        os.remove("coefficients.txt")
    try:
        result = subprocess.run(
            ["sollya", "tmp_interval.sollya"],
            check=True,
            capture_output=True,
            text=True
        )
    except subprocess.CalledProcessError as e:
        return

def print_remez_coeffs(poly):
    print("J1TaylorExtendedSeries {")
    print_double_double("a0: ", poly[0])
    print_double_double("a1: ", poly[1])
    print_double_double("a2: ", poly[2])
    print_double_double("a3: ", poly[3])
    print_double_double("a4: ", poly[4])
    print_double_double("a5: ", poly[5])
    print_double_double("a6: ", poly[6])
    print_double_double("a7: ", poly[7])
    print_double_double("a8: ", poly[8])
    print_double_double("a9: ", poly[9])
    print_double_double("a10: ", poly[10])
    print_double_double("a11: ", poly[11])
    print_double_double("a12: ", poly[12])
    print_double_double("a13: ", poly[13])
    print("c: [")
    for i in range(14, len(poly)):
        coeff = poly[i]
        print(f"{double_to_hex(coeff)},")
    print("],")
    print("},")

degree = 23

print(f"pub(crate) static J1_COEFFS: [J1TaylorExtendedSeries; {len(intervals)}] = [")
for i in range(0, len(intervals)):
    interval = intervals[i]
    call_sollya_on_interval(interval[0], interval[1], interval[2], degree)
    coeffs = load_coefficients(f"coefficients.txt")
    print_remez_coeffs(coeffs)
print("];")