Skip to main content

Y0F_COEFFS

Static Y0F_COEFFS 

Source
pub(crate) static Y0F_COEFFS: [[u64; 18]; 47]
Expand description

Taylor series at different zeros and extremums for Y0 Generated by SageMath:

def compute_intervals(zeros):
    intervals = []
    for i in range(0, len(zeros)):
        if i == 0:
            a = 2 - 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(y0_zeros)
print(intervals)

def build_sollya_script(a, b, zero, deg):
    return f"""
prec = 500;
bessel_y0 = library("./notes/bessel_sollya/cmake-build-release/libbessel_sollya.dylib");
f = bessel_y0(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 [RR(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

degree = 17

print(f"pub(crate) static Y0F_COEFFS: [[u64;{degree + 1}]; {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("[")
    for c in coeffs:
        print(double_to_hex(c) + ",")
    print("],")
print("];")