#!/usr/bin/env python3 """Reproduce the bounded spin-weighted spheroidal-harmonic benchmark. Only (s, ell, m)=(-2,2,2) and the five CASES below are benchmarked. This is not a general branch finder, a quasinormal-mode solver, or an error bound for arbitrary complex c. Every continuation follows the straight segment from c=0, with smaller-step, larger-basis, and independent continued-fraction checks. Decimal parameter strings are passed directly to mpmath. Conventions: exp(-i omega t + i m phi), c=a omega, A(0)=ell(ell+1)-s(s+1), and integral_{-1}^1 |S(x;c)|^2 dx=1. The full-sphere harmonic is S(cos(theta);c) exp(i m phi)/sqrt(2 pi). The real spherical phase is S_j(x;0)=sqrt((2j+1)/32) (1+x)^2 P_{j-2}^{(0,4)}(x). Independent representation: Berti, Cardoso, and Casals, Phys. Rev. D 73, 024013 (2006), corrected arXiv:gr-qc/0511111v4, printed pp. 3-6, Eq. (2.1), Eqs. (2.5)-(2.10), (2.13), and (2.16a)-(2.16b). https://arxiv.org/abs/gr-qc/0511111 Run: python3 spheroidal-harmonics.py --output spheroidal-harmonics-results.json Check an existing result numerically: python3 spheroidal-harmonics.py --check The latter tolerates LAPACK roundoff; it does not overwrite the saved result. All recorded fields are compared, including exact runtime versions and arbitrary-precision decimal strings. A new runtime requires regenerating the result. The check also rejects deliberately corrupted JSON copies in a temporary directory; it never mutates the canonical result. """ from __future__ import annotations import argparse import copy import hashlib import json import math from pathlib import Path import platform import tempfile import mpmath as mp import numpy as np SPIN, ELL, M = -2, 2, 2 CASES = (("0", "0"), ("0.1", "0"), ("0.5", "0"), ("0.3", "0.1"), ("2", "0")) JMAX_VALUES = (8, 12, 16, 24) DEFAULT_OUTPUT = Path(__file__).with_name("spheroidal-harmonics-results.json") def require(condition, message): if not condition: raise RuntimeError(message) def number_pair(value): return {"real": float(np.real(value)), "imag": float(np.imag(value))} def mp_pair(value, digits=65): return {"real": mp.nstr(mp.re(value), digits), "imag": mp.nstr(mp.im(value), digits)} def coupling(j): return np.sqrt((j * j - M * M) * (j * j - SPIN * SPIN) / (j * j * (4 * j * j - 1))) def multiplication_matrix(jmax): js = np.arange(2, jmax + 1) x = np.diag(-SPIN * M / (js * (js + 1))) off_diagonal = np.array([coupling(j) for j in js[1:]]) return x + np.diag(off_diagonal, 1) + np.diag(off_diagonal, -1) def galerkin_matrix(c, jmax): # The extra row retains paths jmax -> jmax+1 -> jmax in P X^2 P. # Squaring the already cropped P X P would omit C_{jmax+1}^2. extended_x = multiplication_matrix(jmax + 1) x = extended_x[:-1, :-1] x_squared = (extended_x @ extended_x)[:-1, :-1] js = np.arange(2, jmax + 1) return np.diag(js * (js + 1) - SPIN * (SPIN + 1)) \ + 2 * c * SPIN * x - c * c * x_squared def matrix_branch(c, jmax=24, max_step=0.05): """Follow the ell=2 branch by overlap, then independently solve its adjoint.""" steps = max(1, int(np.ceil(abs(c) / max_step))) previous = np.zeros(jmax - 1, dtype=complex) previous[0] = 1 minimum_overlap = 1.0 largest_competing_overlap = 0.0 for index in range(1, steps + 1): parameter = c * index / steps matrix = galerkin_matrix(parameter, jmax) if complex(parameter).imag == 0: values, vectors = np.linalg.eigh(matrix.real) else: values, vectors = np.linalg.eig(matrix) overlaps = previous.conj() @ vectors chosen = int(np.argmax(np.abs(overlaps))) candidate = vectors[:, chosen].astype(complex) candidate /= np.linalg.norm(candidate) phase_overlap = np.vdot(previous, candidate) minimum_overlap = min(minimum_overlap, float(abs(phase_overlap))) candidate *= np.exp(-1j * np.angle(phase_overlap)) other = np.delete(np.abs(overlaps), chosen) largest_competing_overlap = max(largest_competing_overlap, float(other.max())) require(abs(phase_overlap) > 0.98 and other.max() < 0.1, "Branch overlap is ambiguous outside this benchmark's safe regime.") previous = candidate value = complex(values[chosen]) # A final b_2 > 0 convention permits comparisons across step sizes. # For complex c this is a smooth unit-norm gauge, not a holomorphic gauge. previous *= np.exp(-1j * np.angle(previous[0])) matrix = galerkin_matrix(c, jmax) left_values, left_vectors = np.linalg.eig(matrix.conj().T) left = left_vectors[:, np.argmin(np.abs(left_values - value.conjugate()))] left /= np.linalg.norm(left) right_residual = float(np.linalg.norm(matrix @ previous - value * previous)) left_residual = float(np.linalg.norm(matrix.conj().T @ left - value.conjugate() * left)) condition = float(1 / abs(np.vdot(left, previous))) # Embedding in two more rows tests omitted couplings of this pentadiagonal H. enlarged = galerkin_matrix(c, jmax + 2) padded = np.pad(previous, (0, 2)) exterior_residual = float(np.linalg.norm(enlarged @ padded - value * padded)) scale = float(np.linalg.norm(matrix, ord=2) + abs(value)) orthogonality = None if complex(c).imag == 0: _, all_vectors = np.linalg.eigh(matrix.real) orthogonality = float(np.linalg.norm(all_vectors.T @ all_vectors - np.eye(jmax - 1), ord=2)) diagnostics = { "jmax": jmax, "steps": steps, "maximum_step": max_step, "A": number_pair(value), "right_residual_l2": right_residual, "left_residual_l2": left_residual, "scaled_right_residual": right_residual / scale, "scaled_left_residual": left_residual / scale, "extended_residual_l2": exterior_residual, "left_right_condition_number": condition, "tail_last_four_l2": float(np.linalg.norm(previous[-4:])), "tail_last_four_max": float(np.max(np.abs(previous[-4:]))), "coefficient_norm_error": float(abs(np.vdot(previous, previous) - 1)), "minimum_step_overlap": minimum_overlap, "largest_competing_step_overlap": largest_competing_overlap, "real_c_basis_orthogonality_error": orthogonality, } return value, previous, diagnostics def recurrence(p, a, c): # k_minus=2, k_plus=0 for the benchmark's s=-2,m=2. alpha = -2 * (p + 1) * (p + 5) beta = p * (p - 1) + 2 * p * (3 - 2 * c) - 6 * c + 4 - c * c - a gamma = 2 * c * p return alpha, beta, gamma def continued_fraction(a, c, cutoff): denominator = recurrence(cutoff, a, c)[1] for p in range(cutoff - 1, 0, -1): alpha, beta, _ = recurrence(p, a, c) gamma_next = recurrence(p + 1, a, c)[2] denominator = beta - alpha * gamma_next / denominator alpha, beta, _ = recurrence(0, a, c) return beta - alpha * recurrence(1, a, c)[2] / denominator def cf_root(c, guess, cutoff, dps): root = mp.findroot(lambda a: continued_fraction(a, c, cutoff), (guess, guess + mp.mpf("0.001")), tol=mp.power(10, -(dps - 10)), maxsteps=50) require(mp.isfinite(root) and abs(continued_fraction(root, c, cutoff)) < mp.power(10, -(dps - 10)), "Continued-fraction root has a nonfinite value or excessive residual.") return root def cf_branch(real, imag, dps=80, cutoff=240, max_step="0.05"): # The initial root is analytic and never comes from a matrix eigenvalue. with mp.workdps(dps): target = mp.mpc(real, imag) steps = max(1, int(mp.ceil(abs(target) / mp.mpf(max_step)))) root = mp.mpc(4) maximum_jump = mp.mpf(0) for index in range(1, steps + 1): c = target * mp.mpf(index) / steps new_root = cf_root(c, root, cutoff, dps) maximum_jump = max(maximum_jump, abs(new_root - root)) require(abs(new_root - root) < mp.mpf("0.5"), "Continued-fraction continuation left the benchmark's safe regime.") root = new_root return root, {"dps": dps, "cutoff": cutoff, "steps": steps, "maximum_step": max_step, "A": mp_pair(root), "absolute_fraction_residual": mp.nstr(abs(continued_fraction(root, target, cutoff)), 8), "largest_step_eigenvalue_change": mp.nstr(maximum_jump, 16)} def jacobi_zero_four(x, highest_degree): """Evaluate Jacobi P_n^(0,4) by the standard three-term recurrence.""" x = np.asarray(x, dtype=float) polynomials = [np.ones_like(x)] if highest_degree > 0: polynomials.append(3 * x - 2) for n in range(2, highest_degree + 1): denominator = 2 * n * (n + 4) * (2 * n + 2) first = (2 * n + 3) * ((2 * n + 4) * (2 * n + 2) * x - 16) second = 2 * (n - 1) * (n + 3) * (2 * n + 4) polynomials.append((first * polynomials[-1] - second * polynomials[-2]) / denominator) return np.array(polynomials) def spherical_basis(x, jmax): js = np.arange(2, jmax + 1) return np.sqrt((2 * js[:, None] + 1) / 32) * (1 + np.asarray(x)) ** 2 \ * jacobi_zero_four(x, jmax - 2) def profile(vector, x): return vector @ spherical_basis(x, len(vector) + 1) def minimal_coefficients(a, c, cutoff=240): """Backward minimal ratios R_p=a_p/a_(p-1), with a_0=1.""" ratios = [mp.mpc(0)] * (cutoff + 2) for p in range(cutoff, 0, -1): alpha, beta, gamma = recurrence(p, a, c) ratios[p] = -gamma / (beta + alpha * ratios[p + 1]) coefficients = [mp.mpc(1)] for p in range(1, cutoff + 1): coefficients.append(coefficients[-1] * ratios[p]) # Bound the discarded finite tail of the actual unnormalized S'' for # S=exp(cx) sum a_p(1+x)^(p+2), uniformly on -1<=x<=1. Product-rule # contributions from the exponential must not be omitted. The bound does # not cover rounding or the uncomputed infinite tail; cutoff/precision # sweeps test those effects separately. tail = mp.mpf(0) keep = len(coefficients) while keep > 2: p = keep - 1 addition = mp.exp(abs(mp.re(c))) * abs(coefficients[p]) * 2 ** p \ * ((p + 2) * (p + 1) + 4 * abs(c) * (p + 2) + 4 * abs(c) ** 2) if tail + addition > mp.power(10, -(mp.mp.dps - 15)): break tail += addition keep -= 1 return coefficients[:keep], tail def series_value_derivatives(coefficients, c, x): y = 1 + x polynomial = mp.mpc(0) derivative = mp.mpc(0) second_derivative = mp.mpc(0) for coefficient in reversed(coefficients): second_derivative = second_derivative * y + 2 * derivative derivative = derivative * y + polynomial polynomial = polynomial * y + coefficient q = y * y * polynomial qp = 2 * y * polynomial + y * y * derivative qpp = 2 * polynomial + 4 * y * derivative + y * y * second_derivative exponential = mp.exp(c * x) return (exponential * q, exponential * (qp + c * q), exponential * (qpp + 2 * c * qp + c * c * q)) def independent_profile_check(real, root, vector): """Compare Jacobi and minimal-series profiles; test the ODE directly.""" with mp.workdps(80): c = mp.mpf(real) coefficients, tail_bound = minimal_coefficients(root, c) raw = lambda x: series_value_derivatives(coefficients, c, x)[0] normalization = mp.sqrt(mp.quad(lambda x: abs(raw(x)) ** 2, [-1, 0, 1])) points = np.linspace(-1, 1, 161) independent = np.array([float(mp.re(raw(mp.mpf(str(x))) / normalization)) for x in points]) matrix_profile = profile(vector, points).real max_difference = float(np.max(np.abs(independent - matrix_profile))) ode_residual = mp.mpf(0) for x_text in ("-0.999", "-0.99", "-0.9", "-0.5", "0", "0.5", "0.9", "0.99", "0.999"): x = mp.mpf(x_text) value, derivative, second = series_value_derivatives(coefficients, c, x) potential = c * c * x * x + 4 * c * x - 2 + root \ - 4 * (1 - x) / (1 + x) terms = ((1 - x * x) * second, -2 * x * derivative, potential * value) scaled = abs(sum(terms)) / (sum(abs(term) for term in terms) + normalization) ode_residual = max(ode_residual, scaled) south, south_derivative, _ = series_value_derivatives(coefficients, c, mp.mpf(-1)) north, north_derivative, _ = series_value_derivatives(coefficients, c, mp.mpf(1)) north_boundary_residual = abs(-2 * north_derivative + (c * c + 4 * c - 2 + root) * north) / normalization near_south = [] for epsilon_text in ("0.0001", "0.000001"): epsilon = mp.mpf(epsilon_text) observed = raw(-1 + epsilon) / epsilon ** 2 limiting = mp.exp(-c) near_south.append({"epsilon": epsilon_text, "relative_error_in_quadratic_pole_limit": mp.nstr(abs(observed / limiting - 1), 14)}) south_errors = [mp.mpf(row["relative_error_in_quadratic_pole_limit"]) for row in near_south] require(all(mp.isfinite(error) and error > 0 for error in south_errors) and south_errors[1] < mp.mpf("1e-5") and south_errors[1] < mp.mpf("0.02") * south_errors[0], "The nontrivial south-pole quadratic limit fails to converge.") require(mp.isfinite(north) and abs(north / normalization) > mp.mpf("1e-3"), "The regular north-pole value is nonfinite or vanishing.") # Convergence of an eigenvalue does not establish convergence of its # eigenfunction. Rebuild and renormalize minimal-series profiles at # independently varied precision and cutoff, including both poles. comparison_points = ("-1", "-0.99", "-0.9", "-0.5", "0", "0.5", "0.9", "0.99", "1") reference_values = [raw(mp.mpf(x)) / normalization for x in comparison_points] profile_sweep = [] for dps, cutoff in ((50, 80), (50, 160), (80, 80), (80, 160), (80, 240)): with mp.workdps(dps): parameter = mp.mpf(real) sweep_root = cf_root(parameter, root + mp.mpf("0.0001"), cutoff, dps) series, discarded_bound = minimal_coefficients(sweep_root, parameter, cutoff) sweep_raw = lambda x: series_value_derivatives(series, parameter, x)[0] sweep_norm = mp.sqrt(mp.quad(lambda x: abs(sweep_raw(x)) ** 2, [-1, 0, 1])) values = [sweep_raw(mp.mpf(x)) / sweep_norm for x in comparison_points] error = max(abs(value - reference) for value, reference in zip(values, reference_values)) require(error < mp.mpf("1e-30"), "Minimal-series profile lacks precision/cutoff convergence.") profile_sweep.append({"dps": dps, "cutoff": cutoff, "series_terms_retained": len(series), "finite_tail_derivative_bound": mp.nstr(discarded_bound, 8), "maximum_normalized_profile_difference_from_80dps_240cutoff": mp.nstr(error, 8)}) require(max_difference < 3e-12, "Independent angular profiles disagree.") require(ode_residual < mp.mpf("1e-50"), "Minimal-series angular ODE residual is too large.") require(north_boundary_residual < mp.mpf("1e-50"), "North-pole regular boundary condition failed.") require(south == 0 and south_derivative == 0, "South-pole quadratic zero failed.") return {"series_terms_retained": len(coefficients), "series_finite_tail_derivative_bound": mp.nstr(tail_bound, 8), "comparison_grid_points": len(points), "matrix_series_profile_max_difference": max_difference, "maximum_scaled_ode_residual": mp.nstr(ode_residual, 8), "north_pole_boundary_residual": mp.nstr(north_boundary_residual, 8), "north_pole_value": mp_pair(north / normalization), "finite_tail_bound_meaning": "Uniform bound on discarded finite terms in the unnormalized S second derivative; excludes rounding and uncomputed infinite tail.", "south_pole_value": mp_pair(south), "south_pole_derivative": mp_pair(south_derivative), "south_pole_approach": near_south, "series_profile_convergence_points": list(comparison_points), "series_profile_convergence_threshold": "1e-30", "series_profile_precision_cutoff_sweep": profile_sweep, "independent_series_norm": mp.nstr(mp.quad(lambda x: abs(raw(x) / normalization) ** 2, [-1, 0, 1]), 20)} def projection_checks(): jmax = 8 extended_x = multiplication_matrix(jmax + 1) x = extended_x[:-1, :-1] difference = (extended_x @ extended_x)[:-1, :-1] - x @ x expected = np.zeros_like(x) expected[-1, -1] = coupling(jmax + 1) ** 2 require(np.allclose(difference, expected, rtol=0, atol=2e-16), "The Galerkin X-squared projection correction is incorrect.") nodes, weights = np.polynomial.legendre.leggauss(160) basis = spherical_basis(nodes, 24) orthogonality = (basis * weights) @ basis.T gram_error = float(np.linalg.norm(orthogonality - np.eye(23), ord=2)) integrated_x_squared = (basis * (weights * nodes ** 2)) @ basis.T ext = multiplication_matrix(25) x_squared_error = float(np.max(np.abs(integrated_x_squared - (ext @ ext)[:-1, :-1]))) require(gram_error < 2e-11 and x_squared_error < 2e-12, "Jacobi quadrature failed an independent basis or X-squared check.") incorrect = galerkin_matrix(2.0, jmax) + 4 * difference correct_value = np.linalg.eigvalsh(galerkin_matrix(2.0, jmax))[0] wrong_value = np.linalg.eigvalsh(incorrect)[0] return {"jmax": jmax, "nonzero_correction_entry": [jmax, jmax], "PX2P_minus_PXP_squared_last_diagonal": float(difference[-1, -1]), "expected_C_jmax_plus_1_squared": float(expected[-1, -1]), "projection_identity_max_error": float(np.max(np.abs(difference - expected))), "jacobi_quadrature_basis_gram_error": gram_error, "jacobi_quadrature_X_squared_max_error": x_squared_error, "c_2_wrong_truncation_eigenvalue_minus_correct": float(wrong_value - correct_value)} def perturbation_checks(): coefficient = float(10 / (9 * np.sqrt(35))) target = np.zeros(23) target[1] = coefficient spherical = np.zeros(23) spherical[0] = 1 rows = [] previous_error = None for c in (0.1, 0.05, 0.025, 0.0125): value, vector, _ = matrix_branch(c) remainder = float(np.linalg.norm(vector.real - spherical - c * target)) derivative_error = remainder / c ratio = None if previous_error is None else previous_error / remainder if ratio is not None: require(3.8 < ratio < 4.2, "First-order profile remainder is not quadratic in c.") eigenvalue_remainder = value.real - (4 - 8 * c / 3) eigenvalue_ratio = None if not rows else rows[-1]["eigenvalue_first_order_remainder"] / eigenvalue_remainder if eigenvalue_ratio is not None: require(3.8 < eigenvalue_ratio < 4.2, "First-order eigenvalue remainder is not quadratic in c.") rows.append({"c": c, "profile_first_order_remainder_l2": remainder, "finite_difference_derivative_error_l2": derivative_error, "previous_remainder_over_current": ratio, "eigenvalue_first_order_remainder": eigenvalue_remainder, "previous_eigenvalue_remainder_over_current": eigenvalue_ratio}) previous_error = remainder return {"A_first_derivative_at_zero": -8 / 3, "S_3_coefficient_in_S_first_derivative": coefficient, "real_c_halving_checks": rows} def run_benchmark(): output = { "schema_version": 2, "source_sha256": hashlib.sha256(Path(__file__).read_bytes()).hexdigest(), "runtime": {"python": platform.python_version(), "numpy": np.__version__, "mpmath": mp.__version__}, "scope": "Only the s=-2, ell=m=2 branch on the five stated straight continuation paths; no QNM or arbitrary-complex-c claim.", "conventions": {"s": SPIN, "ell": ELL, "m": M, "angular_eigenvalue": "A(0)=ell(ell+1)-s(s+1)", "polar_normalization": "integral_-1^1 |S|^2 dx = 1", "full_sphere_harmonic": "S(cos(theta);c) exp(i m phi)/sqrt(2 pi)", "phase": "continuous overlap phase, then b_2 real positive; real S(x;0)=sqrt(5/32)(1+x)^2"}, "primary_source": "https://arxiv.org/abs/gr-qc/0511111 (v4), printed pp. 3-6, Eq. (2.1), Eqs. (2.5)-(2.10), (2.13), (2.16a)-(2.16b)", "projection_checks": projection_checks(), "cases": [], } figure_x = np.linspace(-1, 1, 401) figure_profiles = [] nodes, weights = np.polynomial.legendre.leggauss(160) for real, imag in CASES: c = complex(float(real), float(imag)) matrices = [matrix_branch(c, jmax) for jmax in JMAX_VALUES] value, vector, best = matrices[-1] half_value, half_vector, half_diagnostics = matrix_branch(c, max_step=0.025) root, reference_cf = cf_branch(real, imag) sweep = [] with mp.workdps(80): for dps, cutoff in ((50, 80), (50, 160), (80, 160), (80, 240)): with mp.workdps(dps): # Move off the high-precision answer so findroot must solve # anew rather than return an already-converged initial seed. sweep_root = cf_root(mp.mpc(real, imag), root + mp.mpf("0.0001"), cutoff, dps) residual = abs(continued_fraction(sweep_root, mp.mpc(real, imag), cutoff)) difference = abs(sweep_root - root) require(difference < mp.mpf("1e-38"), "Continued fraction lacks precision/cutoff convergence.") sweep.append({"dps": dps, "cutoff": cutoff, "A": mp_pair(sweep_root, min(65, dps - 5)), "difference_from_80dps_240cutoff": mp.nstr(difference, 8), "absolute_fraction_residual": mp.nstr(residual, 8)}) half_root, half_cf = cf_branch(real, imag, max_step="0.025") cf_step_change = abs(root - half_root) matrix_cf_difference = abs(mp.mpc(value.real, value.imag) - root) require(matrix_cf_difference < mp.mpf("2e-12"), "Matrix and independently continued CF eigenvalues disagree.") require(cf_step_change < mp.mpf("1e-60"), "CF branch changed under step halving.") profile_norm_error = float(abs(np.sum(weights * np.abs(profile(vector, nodes)) ** 2) - 1)) phase_aligned_step_change = float(np.linalg.norm(vector - half_vector)) require(profile_norm_error < 3e-13, "Polar profile normalization failed.") require(abs(value - half_value) < 2e-12 and phase_aligned_step_change < 2e-12, "Matrix continuation changed under step halving.") require(best["right_residual_l2"] < 3e-11 and best["left_residual_l2"] < 3e-11, "A matrix or adjoint eigenvector residual is too large.") require(best["tail_last_four_l2"] < 1e-14, "The final spherical coefficient tail is unresolved.") row = {"c": {"real": real, "imag": imag}, "matrix": best, "matrix_jmax_sweep": [entry[2] for entry in matrices], "matrix_A_change_16_to_24": float(abs(matrices[-2][0] - value)), "continued_fraction": reference_cf, "continued_fraction_sweep": sweep, "matrix_cf_absolute_difference": float(matrix_cf_difference), "matrix_step_halving": {"eigenvalue_change": float(abs(value - half_value)), "phase_aligned_vector_change_l2": phase_aligned_step_change, "steps": half_diagnostics["steps"]}, "cf_step_halving": {"eigenvalue_change": mp.nstr(cf_step_change, 8), "steps": half_cf["steps"]}, "polar_quadrature_norm_error": profile_norm_error, "leading_mixing_coefficients": [{"j": j + 2, "b": number_pair(coefficient)} for j, coefficient in enumerate(vector[:7])], } if imag == "0" and real in ("0.1", "0.5", "2"): row["independent_series_profile_check"] = independent_profile_check(real, root, vector) if imag == "0" and real in ("0", "0.1", "0.5"): figure_profiles.append({"c": real, "S": profile(vector, figure_x).real.tolist()}) output["cases"].append(row) output["perturbation_checks"] = perturbation_checks() predicted_derivative = 10 / (9 * np.sqrt(35)) * spherical_basis(figure_x, 3)[1] difference_quotient = (np.array(figure_profiles[1]["S"]) - np.array(figure_profiles[0]["S"])) / 0.1 output["figure_data"] = {"x": figure_x.tolist(), "profiles": figure_profiles, "difference_quotient_c": "0.1", "difference_quotient": difference_quotient.tolist(), "first_order_derivative": predicted_derivative.tolist(), "all_axes_dimensionless": True} output["checks_passed"] = True validate_diagnostics(output) return output def validate_diagnostics(result): """Apply the same scientific acceptance gates to generated and saved data.""" def bounded(value, upper, name, lower=0): with mp.workdps(80): parsed = mp.mpf(value) require(mp.isfinite(parsed) and mp.mpf(lower) <= parsed < mp.mpf(upper), f"Diagnostic {name} is nonfinite or outside its acceptance interval.") require(result["checks_passed"] is True, "The benchmark did not pass its checks.") for case in result["cases"]: real_parameter = case["c"]["imag"] == "0" for matrix in case["matrix_jmax_sweep"] + [case["matrix"]]: for key in ("right_residual_l2", "left_residual_l2"): bounded(matrix[key], "1e-10", key) for key in ("scaled_right_residual", "scaled_left_residual", "coefficient_norm_error"): bounded(matrix[key], "1e-11", key) bounded(matrix["left_right_condition_number"], "10", "condition number", "0.9999999999") bounded(matrix["minimum_step_overlap"], "1.000000000001", "branch overlap", "0.98") bounded(matrix["largest_competing_step_overlap"], "0.1", "competing overlap") if real_parameter: bounded(matrix["real_c_basis_orthogonality_error"], "1e-11", "real-c orthogonality") bounded(abs(matrix["left_right_condition_number"] - 1), "1e-11", "Hermitian condition number") else: require(matrix["real_c_basis_orthogonality_error"] is None, "Complex-c runs must not claim Hermitian orthogonality.") bounded(case["matrix"]["extended_residual_l2"], "1e-10", "final extended residual") bounded(case["matrix"]["tail_last_four_l2"], "1e-14", "final coefficient tail") bounded(case["matrix_A_change_16_to_24"], "2e-12", "J convergence") bounded(case["matrix_cf_absolute_difference"], "2e-12", "matrix-CF agreement") bounded(case["polar_quadrature_norm_error"], "3e-13", "polar norm") for key in ("eigenvalue_change", "phase_aligned_vector_change_l2"): bounded(case["matrix_step_halving"][key], "2e-12", "matrix step halving") bounded(case["cf_step_halving"]["eigenvalue_change"], "1e-60", "CF step halving") for fraction in [case["continued_fraction"]] + case["continued_fraction_sweep"]: bounded(fraction["absolute_fraction_residual"], f"1e-{fraction['dps'] - 10}", "CF residual") if "difference_from_80dps_240cutoff" in fraction: bounded(fraction["difference_from_80dps_240cutoff"], "1e-38", "CF precision/cutoff convergence") bounded(case["continued_fraction"]["largest_step_eigenvalue_change"], "0.5", "CF branch jump") if "independent_series_profile_check" in case: profile_check = case["independent_series_profile_check"] bounded(profile_check["matrix_series_profile_max_difference"], "3e-12", "matrix-series profile") bounded(profile_check["maximum_scaled_ode_residual"], "1e-50", "angular ODE") bounded(profile_check["north_pole_boundary_residual"], "1e-50", "north boundary") with mp.workdps(80): north = mp.mpc(profile_check["north_pole_value"]["real"], profile_check["north_pole_value"]["imag"]) require(mp.isfinite(north) and abs(north) > mp.mpf("1e-3"), "North-pole value vanishes or is nonfinite.") south_errors = [mp.mpf(row["relative_error_in_quadratic_pole_limit"]) for row in profile_check["south_pole_approach"]] require(len(south_errors) == 2 and all(mp.isfinite(e) and e > 0 for e in south_errors) and south_errors[1] < mp.mpf("1e-5") and south_errors[1] < mp.mpf("0.02") * south_errors[0], "South-pole quadratic-limit convergence failed.") bounded(abs(mp.mpf(profile_check["independent_series_norm"]) - 1), "1e-60", "series norm") for entry in profile_check["series_profile_precision_cutoff_sweep"]: bounded(entry["maximum_normalized_profile_difference_from_80dps_240cutoff"], "1e-30", "profile precision/cutoff convergence") bounded(entry["finite_tail_derivative_bound"], f"1e-{entry['dps'] - 15}", "finite S-second-derivative tail") for entry in result["perturbation_checks"]["real_c_halving_checks"][1:]: for key in ("previous_remainder_over_current", "previous_eigenvalue_remainder_over_current"): bounded(entry[key], "4.2", "perturbative remainder ratio", "3.8") def compare_saved_tree(saved, current, path=()): """Compare every field, never silently ignoring new or missing metadata. Decimal strings are exact because mpmath with the recorded runtime is deterministic. Floating diagnostics have explicit roundoff tolerances; tiny coefficient tails use relative tolerance instead of a blanket absolute tolerance that would erase their significance. """ location = ".".join(map(str, path)) or "root" require(type(saved) is type(current), f"Saved type differs at {location}.") if isinstance(current, dict): require(saved.keys() == current.keys(), f"Saved key set differs at {location}.") for key in current: compare_saved_tree(saved[key], current[key], path + (key,)) elif isinstance(current, list): require(len(saved) == len(current), f"Saved list length differs at {location}.") for index, (old, new) in enumerate(zip(saved, current)): compare_saved_tree(old, new, path + (index,)) elif isinstance(current, float): require(math.isfinite(saved) and math.isfinite(current), f"Nonfinite float at {location}.") field = path[-1] if field == "maximum_step" or (len(path) > 1 and path[-2] == "x"): require(saved == current, f"Saved sampling parameter differs at {location}.") return if len(path) > 1 and path[-2] == "A": # The page prints 12 decimal places; a 2e-12 tolerance could # silently change its table. Keep eigenvalue comparison tighter. rtol, atol = 0, 1e-13 elif str(field).startswith("tail_last_four_"): rtol, atol = 1e-3, 1e-60 elif field == "left_right_condition_number": rtol, atol = 1e-12, 1e-12 elif any(word in str(field) for word in ("residual", "error", "difference", "change")): # Diagnostics are reported measurements, sometimes far below # 1e-13. Preserve their scale instead of accepting an absolute # tolerance larger than the reported value itself. Different # hardware may require regenerating its own measured report. rtol, atol = 1e-3, 1e-30 else: rtol, atol = 1e-11, 2e-12 require(abs(saved - current) <= atol + rtol * abs(current), f"Saved numerical value differs at {location}.") else: require(saved == current, f"Saved metadata or decimal value differs at {location}.") def check_saved(saved, current): compare_saved_tree(saved, current) validate_diagnostics(saved) validate_diagnostics(current) def saved_result_poison_checks(current): """Prove corrupt saved reports fail, using only isolated temporary JSONs.""" mutations = [ ("right-residual", ("cases", 1, "matrix", "right_residual_l2"), 1.0), ("realistic-right-residual", ("cases", 3, "matrix", "right_residual_l2"), 9e-14), ("realistic-profile-disagreement", ("cases", 1, "independent_series_profile_check", "matrix_series_profile_max_difference"), 9e-14), ("displayed-eigenvalue", ("cases", 1, "matrix", "A", "real"), current["cases"][1]["matrix"]["A"]["real"] + 1e-12), ("condition-number", ("cases", 3, "matrix", "left_right_condition_number"), 1000000.0), ("polar-norm", ("cases", 1, "polar_quadrature_norm_error"), 0.2), ("CF-sweep", ("cases", 3, "continued_fraction_sweep", 0, "difference_from_80dps_240cutoff"), "1.0"), ("profile-sweep", ("cases", 1, "independent_series_profile_check", "series_profile_precision_cutoff_sweep", 0, "maximum_normalized_profile_difference_from_80dps_240cutoff"), "1.0"), ("runtime", ("runtime", "numpy"), "unverified-runtime"), ("extended-residual", ("cases", 4, "matrix", "extended_residual_l2"), 0.1), ("J-convergence", ("cases", 4, "matrix_A_change_16_to_24"), 0.1), ("orthogonality", ("cases", 1, "matrix", "real_c_basis_orthogonality_error"), 0.1), ("south-limit", ("cases", 1, "independent_series_profile_check", "south_pole_approach", 1, "relative_error_in_quadratic_pole_limit"), "1.0"), ("tiny-tail", ("cases", 1, "matrix", "tail_last_four_l2"), 1e-20), ("convention", ("conventions", "polar_normalization"), "unnormalized"), ] rejected = [] with tempfile.TemporaryDirectory(prefix="bh-spheroidal-poison-") as directory: for name, path, value in mutations: poisoned = copy.deepcopy(current) target = poisoned for key in path[:-1]: target = target[key] target[path[-1]] = value file = Path(directory) / f"{name}.json" file.write_text(json.dumps(poisoned, allow_nan=False)) try: check_saved(json.loads(file.read_text()), current) except RuntimeError: rejected.append(name) else: raise RuntimeError(f"Saved-result check accepted poison case {name}.") # Structural omissions must fail even when all retained values agree. incomplete = copy.deepcopy(current) del incomplete["cases"][0]["matrix"]["left_residual_l2"] file = Path(directory) / "missing-field.json" file.write_text(json.dumps(incomplete, allow_nan=False)) try: check_saved(json.loads(file.read_text()), current) except RuntimeError: rejected.append("missing-field") else: raise RuntimeError("Saved-result check accepted a missing field.") return rejected def main(): parser = argparse.ArgumentParser(description=__doc__) parser.add_argument("--output", type=Path, default=DEFAULT_OUTPUT, help="Result JSON path; default is beside this script.") parser.add_argument("--check", action="store_true", help="Recompute and compare the existing --output without overwriting it.") args = parser.parse_args() result = run_benchmark() if args.check: check_saved(json.loads(args.output.read_text()), result) rejected = saved_result_poison_checks(result) print(f"Spheroidal benchmark: all scientific checks, complete saved-result comparison, and {len(rejected)} temporary-JSON poison checks passed.") print("Rejected poison cases: " + ", ".join(rejected)) else: args.output.parent.mkdir(parents=True, exist_ok=True) args.output.write_text(json.dumps(result, indent=2, allow_nan=False) + "\n") print(f"Spheroidal benchmark: all scientific checks passed; wrote {args.output}.") if __name__ == "__main__": main()