pub(crate) static J0F_COEFFS: [[u64; 14]; 47]Expand description
Taylor coefficients for J0 [0; 74.62]
Generated by SageMath with 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(j0_zeros)
# print(intervals)
def build_sollya_script(a, b, zero, deg):
return f"""
prec = 500;
bessel_j0 = library("./notes/bessel_sollya/cmake-build-release/libbessel_sollya.dylib");
f = bessel_j0(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 = 13
print(f"pub(crate) static J0F_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("];")