back to table · edit · history · where entries came from · files · download
49201 bytes, as of the version from 2026-11-02 11:01 (current). Recorded here, not run.
"""Offline independent Hessian checks; see README.md for inputs and limits.
Ordinary values use Cardano roots of the ORIGINAL cubic and Jensen quadrature.
Modified values use the corrected hypergeometric identity, never main periods.
Numerical checks support agreement, not certified analytic error bounds.
"""
import argparse
from contextlib import contextmanager
from datetime import datetime, timezone
from fractions import Fraction
import hashlib
import inspect
from itertools import permutations
import json
import math
import os
from pathlib import Path
import platform
import re
import signal
import sys
import time
import mpmath as mp
HERE = Path(__file__).resolve().parent
ROOT = HERE.parent
REQUIRED_DIGITS = 110
EVALUATOR_SHA256 = hashlib.sha256(Path(__file__).read_bytes()).hexdigest()
DECIMAL = re.compile(r"[+-]?(?:\d+(?:\.\d*)?|\.\d+)(?:[eE][+-]?\d+)?\Z")
RATIONAL = re.compile(r"[+-]?\d+(?:/[1-9]\d*)?\Z")
SIGNS = {
0: ((1, 1, 1, 1), (1, 1, 1, 1), (1, 1, 1, 1), (1, 1, 1, 1)),
1: ((1, 1, 1, 1), (1, 1, 1, 1), (1, 1, -1, -1), (1, 1, -1, -1)),
2: ((1, 1, 1, 1), (1, -1, -1, 1), (1, -1, 1, -1), (1, 1, -1, -1)),
3: ((1, 1, 1, 1), (1, -1, 1, -1), (1, 1, 1, 1), (1, -1, 1, -1)),
4: ((1, 1, 1, 1), (1, -1, 1, -1), (1, 1, -1, -1), (1, -1, -1, 1)),
5: ((1, 1, 1, 1), (1, -1, -1, 1), (1, -1, -1, 1), (1, 1, 1, 1)),
6: ((1, 1, 1, 1), (1, 1, -1, -1), (1, -1, 1, -1), (1, -1, -1, 1)),
}
EXPANSION_SIGNS = {
0: (-1, -1, -1, -1, -1, -1, 1),
1: (-1, 1, 1, 1, 1, -1, -1),
2: (1, -1, 1, 1, -1, 1, 1),
3: (1, -1, 1, 1, -1, 1, -1),
4: (1, 1, -1, -1, 1, 1, 1),
5: (1, 1, -1, -1, 1, 1, -1),
6: (-1, -1, -1, -1, -1, -1, -1),
}
class CheckFailure(ValueError):
pass
def require(condition, message):
if not condition:
raise CheckFailure(message)
def sha256(path):
return hashlib.sha256(Path(path).read_bytes()).hexdigest()
def object_hash(value):
return hashlib.sha256(json.dumps(value, sort_keys=True).encode()).hexdigest()
def save_json(path, value):
path = Path(path).resolve()
require(path.is_relative_to(HERE), "Output must remain inside independent/")
path.parent.mkdir(parents=True, exist_ok=True)
temporary = path.with_suffix(path.suffix + ".tmp")
temporary.write_text(json.dumps(value, indent=2, allow_nan=False) + "\n")
temporary.replace(path)
def rational(value):
require(type(value) in (str, int), "Exact rationals must be strings or integers")
require(bool(RATIONAL.fullmatch(str(value))), "Not an exact rational: " + repr(value))
return Fraction(value)
def decimal(value):
require(isinstance(value, str) and bool(DECIMAL.fullmatch(value)),
"A finite decimal STRING is required: " + repr(value))
result = mp.mpf(value)
require(mp.isfinite(result), "Nonfinite decimal")
return result
def trim(polynomial):
polynomial = list(polynomial)
while len(polynomial) > 1 and not polynomial[-1]:
polynomial.pop()
return polynomial
def polynomial_value(polynomial, argument):
result = Fraction(0)
for coefficient in reversed(polynomial):
result = result * argument + coefficient
return result
def remainder(dividend, divisor):
dividend, divisor = trim(dividend), trim(divisor)
while dividend != [0] and len(dividend) >= len(divisor):
shift = len(dividend) - len(divisor)
multiplier = dividend[-1] / divisor[-1]
for index, coefficient in enumerate(divisor):
dividend[index + shift] -= multiplier * coefficient
dividend = trim(dividend)
return dividend
def sturm_sequence(polynomial):
polynomial = [Fraction(coefficient) for coefficient in polynomial]
derivative = [index * coefficient for index, coefficient in enumerate(polynomial)][1:]
sequence = [polynomial, derivative]
while sequence[-1] != [0]:
residual = [-coefficient for coefficient in remainder(sequence[-2], sequence[-1])]
if residual == [0]:
break
scale = abs(residual[-1])
sequence.append([coefficient / scale for coefficient in residual])
require(len(sequence[-1]) == 1, "Argument polynomial is not square-free")
return sequence
def variations(sequence, argument):
signs = []
for polynomial in sequence:
value = polynomial_value(polynomial, argument)
if value:
signs.append(1 if value > 0 else -1)
return sum(left != right for left, right in zip(signs, signs[1:]))
def validate_argument(item):
require(isinstance(item, dict), "Argument must be an object")
require(isinstance(item.get("key"), str) and item["key"], "Missing argument key")
polynomial = item.get("polynomial")
require(isinstance(polynomial, list) and 2 <= len(polynomial) <= 5,
"Polynomial must have degree 1 through 4")
require(all(type(coefficient) is int for coefficient in polynomial),
"Polynomial coefficients must be exact JSON integers")
require(polynomial[-1] != 0, "Polynomial has a leading zero")
require(type(item.get("degree")) is int and item["degree"] == len(polynomial) - 1,
"Declared degree disagrees with polynomial")
require(type(item.get("modified")) is bool, "modified must be a boolean")
interval = item.get("isolating_interval")
require(isinstance(interval, list) and len(interval) == 2, "Missing isolating interval")
lower, upper = map(rational, interval)
require(lower <= upper, "Reversed isolating interval")
sequence = sturm_sequence(polynomial)
if lower == upper:
require(item["degree"] == 1 and polynomial_value(polynomial, lower) == 0,
"Collapsed interval is permitted only for an exact linear root")
else:
require(polynomial_value(polynomial, lower) != 0 and
polynomial_value(polynomial, upper) != 0, "Root at isolating endpoint")
require(variations(sequence, lower) - variations(sequence, upper) == 1,
"Interval does not isolate exactly one real root")
return lower, upper
def resolve_argument(item, working):
lower, upper = validate_argument(item)
polynomial = item["polynomial"]
if item["degree"] == 1:
lower = upper = -Fraction(polynomial[0], polynomial[1])
target = Fraction(1, 10 ** (working + 25))
iterations = 0
while lower != upper:
scale = min(abs(lower), abs(upper))
if lower * upper > 0 and upper - lower < scale * target:
break
midpoint = (lower + upper) / 2
middle_value = polynomial_value(polynomial, midpoint)
if middle_value == 0:
require(item["degree"] == 1, "Nonminimal polynomial has a rational root")
lower = upper = midpoint
break
if polynomial_value(polynomial, lower) * middle_value < 0:
upper = midpoint
else:
lower = midpoint
iterations += 1
require(iterations <= 10000, "Exact root refinement exceeded safety limit")
midpoint = (lower + upper) / 2
inside = -1 < midpoint < 27
require(inside == item["modified"], "modified flag disagrees with exact root location")
require(lower == upper or not any(lower <= boundary <= upper for boundary in (-1, 0, 27)),
"Root enclosure does not resolve a branch boundary")
with mp.workdps(working + 20):
value = mp.mpf(midpoint.numerator) / midpoint.denominator
evidence = {
"lower": str(lower), "upper": str(upper), "bisections": iterations,
"exact_rational": str(lower) if lower == upper else None,
"method": "rational Sturm isolation and rational bisection",
}
return value, evidence
def agreement(actual, expected, label, required=REQUIRED_DIGITS):
require(mp.isfinite(actual) and mp.isfinite(expected), label + ": nonfinite value")
require(actual != 0 and expected != 0, label + ": nonzero values required")
require(mp.sign(actual) == mp.sign(expected), label + ": sign mismatch")
relative = abs(actual - expected) / max(abs(actual), abs(expected))
require(relative < mp.power(10, -required),
label + ": insufficient significant-digit agreement (relative error " +
mp.nstr(relative, 12) + ")")
return {
"relative_error": mp.nstr(relative, 15),
"agreement_digits": mp.nstr(-mp.log10(relative), 12) if relative else "equal at working precision",
"required_digits_strictly_greater_than": required,
}
def check_stored(stored, actual, exact=False):
if exact:
require(type(stored) is int and stored == 0 or stored == "0",
"The exact modified zero must not be stored as a rounded decimal")
require(actual == 0, "Exact-zero mismatch")
return {"exact_zero": True}
require(isinstance(stored, str) and ("." in stored or "e" in stored.lower()),
"A nonzero stored approximation must be a decimal string")
centre = decimal(stored)
mantissa, *exponent = re.split("[eE]", stored)
exponent = int(exponent[0]) if exponent else 0
places = len(mantissa.partition(".")[2])
unit = mp.power(10, exponent - places)
distance = abs(actual - centre) / unit
require(distance <= 1, "Independent value lies outside the stored one-ULP interval")
return {"distance_in_stored_ulps": mp.nstr(distance, 15)}
def cubic_jensen(argument, working):
with mp.workdps(working + 40):
mp.mp._tanh_sinh.clear()
parameter = mp.mpf(argument)
cube_root = mp.sign(parameter) * mp.root(abs(parameter), 3)
gap = 3 - cube_root
near_boundary = 0 < gap < mp.mpf("0.01")
root_scale = mp.sqrt(cube_root / 3) if near_boundary else None
unity = mp.exp(2 * mp.pi * mp.j / 3)
calls = 0
max_residual = mp.mpf(0)
def integrand(angle):
nonlocal calls, max_residual
calls += 1
coordinate = mp.exp(mp.j * angle)
linear = -cube_root * coordinate
constant = coordinate ** 3 + 1
if parameter == 0:
return mp.log(max(mp.mpf(1), abs(constant)))
if parameter == 27:
roots = (-coordinate - 1, -unity * coordinate - unity ** 2,
-unity ** 2 * coordinate - unity)
elif near_boundary:
ratio = mp.cos(3 * angle / 2) / root_scale ** 3
phase = mp.exp(mp.j * angle / 2)
if ratio >= 1:
auxiliary = mp.acosh(ratio) / 3
real_part = root_scale * mp.cosh(auxiliary)
imaginary_part = root_scale * mp.sqrt(3) * mp.sinh(auxiliary)
real_roots = (-2 * real_part, real_part + mp.j * imaginary_part,
real_part - mp.j * imaginary_part)
else:
auxiliary = mp.acos(ratio) / 3
real_part = root_scale * mp.cos(auxiliary)
offset = root_scale * mp.sqrt(3) * mp.sin(auxiliary)
real_roots = (-2 * real_part, real_part + offset, real_part - offset)
roots = tuple(phase * root for root in real_roots)
else:
discriminant = mp.sqrt(constant ** 2 / 4 + linear ** 3 / 27)
choices = (-constant / 2 + discriminant, -constant / 2 - discriminant)
first = mp.root(max(choices, key=abs), 3)
require(first != 0, "Cardano divided by zero")
second = -linear / (3 * first)
roots = (first + second, unity * first + unity ** 2 * second,
unity ** 2 * first + unity * second)
if calls <= 12 or calls % 127 == 0:
for root in roots:
scale = abs(root) ** 3 + abs(linear * root) + abs(constant)
residual = abs(root ** 3 + linear * root + constant) / max(1, scale)
max_residual = max(max_residual, residual)
require(residual < mp.power(10, -working + 12), "Cardano root residual too large")
return mp.fsum(mp.log(max(mp.mpf(1), abs(root))) for root in roots)
endpoint = mp.pi / 3
points = [mp.mpf(0), endpoint]
if -1 < parameter < 27:
crossing = mp.acos((cube_root - 1) / 2)
points.extend((crossing / 3, min(2 * crossing / 3, (2 * mp.pi - 2 * crossing) / 3)))
if 0 < parameter < 27:
points.append(2 * mp.acos(mp.sqrt(parameter / 27)) / 3)
geometric_points = 0
if near_boundary:
centers = sorted(set(points))
near_end = min(endpoint, 4 * max(point for point in centers if point < endpoint))
step = min(upper - lower for lower, upper in zip(centers, centers[1:])) / 4
while step < near_end:
for center in centers:
if center < near_end:
for candidate in (center - step, center + step):
if 0 < candidate < near_end:
points.append(candidate)
geometric_points += 1
step *= 4
step = near_end
while step < endpoint:
points.append(step)
geometric_points += 1
step *= 4
points = sorted(set(points))
unique = [points[0]]
for point in points[1:]:
if point - unique[-1] > mp.eps * 100:
unique.append(point)
split = [unique[0]]
for lower, upper in zip(unique, unique[1:]):
split.extend(((lower + upper) / 2, upper))
tolerance = mp.power(10, -working + 10)
def integrate(lower, upper, depth=0):
value, error = mp.quad(integrand, [lower, upper], error=True, maxdegree=8)
require(mp.isfinite(value) and mp.isfinite(error), "Nonfinite Jensen quadrature")
if error < tolerance * (upper - lower) / endpoint:
return value, error
require(depth < 24, "Cubic Jensen quadrature did not converge: " + json.dumps({
"lower": mp.nstr(lower, 30), "upper": mp.nstr(upper, 30),
"width": mp.nstr(upper - lower, 12), "estimate": mp.nstr(error, 12),
"target": mp.nstr(tolerance * (upper - lower) / endpoint, 12),
"gap_3_minus_k": mp.nstr(gap, 15)}))
middle = (lower + upper) / 2
left, left_error = integrate(lower, middle, depth + 1)
right, right_error = integrate(middle, upper, depth + 1)
return left + right, left_error + right_error
pieces = [integrate(lower, upper) for lower, upper in zip(split, split[1:])]
value = 3 * mp.fsum(piece[0] for piece in pieces) / mp.pi
error = 3 * mp.fsum(piece[1] for piece in pieces) / mp.pi
require(mp.isfinite(value) and value > 0, "Ordinary Mahler measure must be finite and positive")
return value, {
"method": "original Hessian cubic Cardano-root Jensen quadrature",
"extra_working_digits": 40,
"arithmetic_working_digits": working + 40,
"geometric_breakpoints": geometric_points,
"near_boundary_root_coordinates": "y=exp(i*theta/2)*u; u^3-k*u+2*cos(3*theta/2)=0" if near_boundary else None,
"calls": calls, "quadrature_error_estimate": mp.nstr(error, 12),
"maximum_sampled_cubic_relative_residual": mp.nstr(max_residual, 12),
"error_estimate_is_rigorous": False,
}
def modified_boundary_series(argument, working):
"""Integrate the zero-balanced 2F1 log series from the t=27 boundary."""
with mp.workdps(working + 20):
parameter = mp.mpf(argument)
complement = (27 - parameter) / 27
require(0 < complement <= mp.mpf(5) / 27, "Boundary series requires 22 <= t < 27")
logarithm = mp.log(complement)
branches = []
for first, second in ((mp.mpf(1) / 3, mp.mpf(2) / 3), (mp.mpf(2) / 3, mp.mpf(1) / 3)):
branches.append({"first": first, "second": second, "hyper": [mp.mpf(1)],
"binomial": [mp.mpf(1)],
"constant": [2 * mp.digamma(1) - 2 * mp.digamma(first)]})
total = mp.mpf(0)
power = complement
target = mp.power(10, -working - 10)
for index in range(2000):
constant_coefficient = mp.mpf(0)
logarithmic_coefficient = mp.mpf(0)
for branch in branches:
if index:
branch["hyper"].append(branch["hyper"][-1] * ((branch["first"] + index - 1) / index) ** 2)
branch["binomial"].append(branch["binomial"][-1] * (branch["second"] + index - 1) / index)
previous = branch["constant"][-1] / branch["hyper"][-2]
difference = previous + 2 / mp.mpf(index) - 2 / (branch["first"] + index - 1)
branch["constant"].append(branch["hyper"][-1] * difference)
constant_coefficient += mp.fsum(branch["binomial"][index - offset] * branch["constant"][offset]
for offset in range(index + 1))
logarithmic_coefficient += mp.fsum(branch["binomial"][index - offset] * branch["hyper"][offset]
for offset in range(index + 1))
term = power / (index + 1) * (constant_coefficient + logarithmic_coefficient / (index + 1)
- logarithmic_coefficient * logarithm)
total += term
if index >= 8 and abs(term) / (1 - complement) < target:
value = -6 * dirichlet_control() + total / (2 * mp.pi * mp.sqrt(3))
return value, {"method": "integrated zero-balanced 2F1 Frobenius log series at t=27",
"extra_working_digits": 20, "arithmetic_working_digits": working + 20,
"terms": index + 1, "last_term": mp.nstr(term, 15),
"boundary": "Ntilde(27-)=-6*mu3", "source": "DLMF 15.8.10 (m=0); Samart Theorem 1",
"error_estimate_is_rigorous": False}
power *= complement
raise CheckFailure("Modified boundary log series exceeded term limit")
def modified_hypergeometric(argument, working):
if 22 <= argument < 27:
return modified_boundary_series(argument, working)
with mp.workdps(working):
parameter = mp.mpf(argument)
require(-1 < parameter < 27, "Modified period requires -1 < t < 27")
if parameter == 0:
return mp.mpf(0), {
"method": "exact continuous value of the corrected convergent identity",
"exact": True,
"arithmetic_working_digits": working,
"proof": "Both 3F2 series are analytic at t=0 and equal 1 there; their prefactors are k and k^2. Continuity of I-2J gives the exact value zero.",
}
cube_root = mp.sign(parameter) * mp.root(abs(parameter), 3)
third = mp.mpf(1) / 3
first_gamma = mp.gamma(third) ** 3
second_gamma = mp.gamma(2 * third) ** 3
if parameter < 0:
extra = max(30, int(-mp.log10(abs(parameter))) + 30)
with mp.extradps(extra):
third = mp.mpf(1) / 3
refined_cube_root = -mp.root(-parameter, 3)
value = mp.re(mp.log(mp.mpc(refined_cube_root)) - 2 / parameter *
mp.hyper([4 * third, 5 * third, 1, 1], [2, 2, 2], 27 / parameter))
method = "Theorem 1 real part of 4F3 continuation; avoids inconsistent printed negative-argument 3F2 prefactor"
error = None
elif parameter < 22:
first = mp.hyper([third] * 3, [2 * third, 4 * third], parameter / 27)
second = mp.hyper([2 * third] * 3, [4 * third, 5 * third], parameter / 27)
value = -(first_gamma * cube_root * first +
second_gamma * cube_root ** 2 * second / 2) / (4 * mp.pi ** 2)
method = "corrected convergent pair of 3F2 series"
error = None
else:
def derivative(coordinate):
argument = coordinate ** 3 / 27
return -(first_gamma * mp.hyp2f1(third, third, 2 * third, argument) +
second_gamma * coordinate * mp.hyp2f1(2 * third, 2 * third, 4 * third, argument)) / (4 * mp.pi ** 2)
points = [mp.mpf(0)]
exponent = 1
while 3 - mp.power(2, -exponent) < cube_root:
points.append(3 - mp.power(2, -exponent))
exponent += 1
require(exponent < 10000, "Hypergeometric path too close to the singularity")
points.append(cube_root)
value, error = mp.quad(derivative, points, method="gauss-legendre", maxdegree=8, error=True)
require(mp.isfinite(error) and error < mp.power(10, -working + 12) * abs(value),
"Integrated hypergeometric derivative did not converge")
method = "integral of the differentiated corrected 3F2 series (two 2F1 functions)"
require(mp.isfinite(value) and value != 0, "Nonfinite or spurious zero modified period")
require(mp.sign(value) == -mp.sign(parameter), "Modified-period sign convention failed")
return value, {"method": method,
"extra_working_digits": extra if parameter < 0 else 0,
"arithmetic_working_digits": working + extra if parameter < 0 else working,
"source_caveat": "Printed post-Theorem-1 3F2 coefficient disagrees with Theorem 1 by a factor of four for negative k" if parameter < 0 else None,
"quadrature_error_estimate": mp.nstr(error, 12) if error is not None else None,
"error_estimate_is_rigorous": False}
def dirichlet_control():
third = mp.mpf(1) / 3
lvalue = (mp.zeta(2, third) - mp.zeta(2, 2 * third)) / 9
return 3 * mp.sqrt(3) * lvalue / (4 * mp.pi)
def exterior_series(parameter):
require(abs(parameter) > 27, "Exterior constant-term series needs |t|>27")
ratio_bound = 27 / abs(parameter)
term = 2 / parameter
total = term
index = 1
while True:
term *= mp.mpf(3 * index * (3 * index + 1) * (3 * index + 2)) / ((index + 1) ** 3 * parameter)
index += 1
total += term
tail = abs(term) * ratio_bound / (1 - ratio_bound)
if tail < mp.power(10, -mp.mp.dps + 6):
return mp.log(abs(parameter)) / 3 - total
require(index < 1000000, "Exterior series exceeded term limit")
def get_pari():
import numberdb.sage
from sage.libs.pari import pari
return pari
def elliptic_control(parameter, kind, working):
pari = get_pari()
pari.set_real_precision(working + 15)
precision = math.ceil((working + 15) * math.log2(10))
with mp.workdps(working):
zeta_check = decimal(str(pari.lfun(1, 2, precision=precision)))
control = agreement(zeta_check, mp.pi ** 2 / 6, "PARI zeta(2) nonzero control")
if kind == "modified" and parameter == 8:
invariants, conductor, multiplier = [2, 0, 1, 0, 0], 19, Fraction(-3)
source = "Samart 2301.05390v3, Theorem 2, E_2: y^2+2xy+y=x^3 (19a3)"
else:
catalogue = {
-216: (27, Fraction(3)), -27: (54, Fraction(1)),
-8: (35, Fraction(1)), -1: (14, Fraction(2)),
32: (20, Fraction(8, 3)), 54: (36, Fraction(3, 2)),
125: (14, Fraction(7)),
}
require(kind == "ordinary" and parameter in catalogue, "No proven elliptic control specified")
conductor, multiplier = catalogue[parameter]
invariants = [0, -27 * parameter ** 2, 0, 216 * parameter ** 3 * (parameter - 27),
-432 * parameter ** 4 * (parameter - 27) ** 2]
source = "Samart 2301.05390v3, Table 1; Samart-Tao (1.7) rational C_t model"
curve = pari.ellinit(invariants)
require(int(curve.ellglobalred()[0]) == conductor, "Unexpected elliptic conductor")
root_number = int(curve.ellrootno())
require(root_number == 1, "The cited rank-zero control has an unexpected root number")
lvalue = decimal(str(pari.lfun(curve, 2, precision=precision)))
require(lvalue > 0, "Nonpositive elliptic L(E,2)")
derivative = root_number * conductor * lvalue / (4 * mp.pi ** 2)
value = mp.mpf(multiplier.numerator) / multiplier.denominator * derivative
return value, {
"source": source, "identity_status": "proved in cited literature",
"pari_version": str(pari("version()")),
"ainvariants": invariants, "conductor": conductor, "root_number": root_number,
"L2": mp.nstr(lvalue, working), "multiplier_of_Lprime0": str(multiplier),
"functional_equation": "L'(E,0)=root_number*N*L(E,2)/(4*pi^2)",
"pari_zeta2_control": control,
}
def signed_determinant(values, coefficients, variant):
require(type(variant) is int and variant in SIGNS, "Unknown determinant sign variant")
require(len(values) == len(coefficients) == 4, "Quartic determinant requires four entries")
require(all(type(coefficient) is int and coefficient > 0 for coefficient in coefficients),
"Determinant coefficients must be positive integers")
weighted = [mp.sqrt(coefficient) * value for coefficient, value in zip(coefficients, values)]
matrix = mp.matrix([[SIGNS[variant][row][column] * weighted[row ^ column]
for column in range(4)] for row in range(4)])
eliminated = mp.det(matrix)
expanded = mp.fsum(
(-1 if sum(order[left] > order[right] for left in range(4) for right in range(left + 1, 4)) % 2 else 1)
* mp.fprod(matrix[row, column] for row, column in enumerate(order))
for order in permutations(range(4))
)
check = agreement(eliminated, expanded, "LU determinant versus 24-term Leibniz expansion")
pairs = ((0, 1), (0, 2), (0, 3), (1, 2), (1, 3), (2, 3))
signs = EXPANSION_SIGNS[variant]
algebraic = mp.fsum(value ** 4 for value in weighted)
algebraic += 2 * mp.fsum(sign * weighted[first] ** 2 * weighted[second] ** 2
for sign, (first, second) in zip(signs[:6], pairs))
algebraic += 8 * signs[6] * mp.fprod(weighted)
algebraic_check = agreement(eliminated, algebraic, "LU determinant versus explicit quartic polynomial")
return eliminated, {"matrix": [[mp.nstr(matrix[row, column], mp.mp.dps) for column in range(4)]
for row in range(4)], "expansion_check": check,
"algebraic_polynomial_check": algebraic_check,
"variant": variant, "absolute_value_taken": False}
def g4_jensen(working):
"""Independent Jensen integral of (1+x)(1+y)(x+y)-4xy."""
with mp.workdps(working):
def integrand(angle):
coordinate = mp.exp(mp.j * angle)
leading = 1 + coordinate
linear = (1 + coordinate) ** 2 - 4 * coordinate
constant = coordinate * leading
discriminant = mp.sqrt(linear ** 2 - 4 * leading * constant)
roots = ((-linear + discriminant) / (2 * leading),
(-linear - discriminant) / (2 * leading))
return mp.log(abs(leading)) + mp.fsum(mp.log(max(1, abs(root))) for root in roots)
crossing = 2 * mp.acos((mp.sqrt(5) - 1) / 2)
points = [mp.mpf(0), crossing / 2, crossing, (crossing + mp.pi) / 2, mp.pi]
value, error = mp.quad(integrand, points, error=True, method="tanh-sinh", maxdegree=9)
require(mp.isfinite(value) and mp.isfinite(error) and error < mp.power(10, -working + 12),
"Independent g(4) quadrature did not converge")
return value / mp.pi, {"method": "direct quadratic-root Jensen integral of the different Rogers-Zudilin polynomial",
"quadrature_error_estimate": mp.nstr(error / mp.pi, 12)}
def provenance(catalogue=None, input_path=None):
paths = [Path(__file__), HERE / "test_verify.py", ROOT / "sources/samart-tao.html",
ROOT / "sources/data.pdf", HERE / "sources/samart.html", HERE / "sources/rogers-zudilin.html"]
if catalogue:
paths.append(Path(catalogue))
if input_path and Path(input_path).is_file():
paths.append(Path(input_path))
try:
from sage.version import version as sage_version
except ImportError:
sage_version = None
return {
"utc": datetime.now(timezone.utc).isoformat(), "python": platform.python_version(),
"sage": sage_version,
"loaded_evaluator_sha256": EVALUATOR_SHA256,
"mpmath": mp.__version__, "hashes": {str(path): sha256(path) for path in paths if path.exists()},
"docker_image_id": os.environ.get("INDEPENDENT_DOCKER_IMAGE_ID"),
"rigour": "heuristic (agreement-checked); argument isolation and the modified zero are exact",
"independence": "No import of shared period engine or shared generators; no network/API calls",
}
@contextmanager
def deadline(seconds):
def expired(signum, frame):
raise TimeoutError("Independent item exceeded time limit; no passing result recorded")
previous = signal.signal(signal.SIGALRM, expired)
signal.setitimer(signal.ITIMER_REAL, seconds)
try:
yield
finally:
signal.setitimer(signal.ITIMER_REAL, 0)
signal.signal(signal.SIGALRM, previous)
def load_catalogue(path):
catalogue = json.loads(Path(path).read_text())
arguments = {}
for item in catalogue["arguments"]:
validate_argument(item)
require(item["key"] not in arguments, "Duplicate catalogue key")
arguments[item["key"]] = item
orbits = {}
for orbit in catalogue.get("quartic_orbits", []):
row = str(orbit["source_row"])
require(row not in orbits, "Duplicate orbit source_row")
require(len(orbit["keys"]) == 4 and len(set(orbit["keys"])) == 4, "Repeated/missing quartic keys")
entries = [arguments[key] for key in orbit["keys"]]
require(all(entry["degree"] == 4 for entry in entries), "Nonquartic orbit argument")
require(all(entry["polynomial"] == entries[0]["polynomial"] for entry in entries),
"Orbit arguments have different polynomials")
require(type(orbit["variant"]) is int and orbit["variant"] in SIGNS, "Invalid matrix variant")
require(len(orbit["coefficients"]) == 4 and
all(type(value) is int and value > 0 for value in orbit["coefficients"]), "Invalid matrix coefficients")
require(rational(orbit["r"]) != 0 and rational(orbit["s"]) != 0, "Source L-value factors must be nonzero")
refined = [resolve_argument(entry, 40)[1] for entry in entries]
require(all(Fraction(left["upper"]) < Fraction(right["lower"])
for left, right in zip(refined, refined[1:])), "Orbit keys are not ordered t1<t2<t3<t4")
orbits[row] = orbit
return arguments, orbits
def read_main_result(directory, item):
filename = hashlib.sha256(item["key"].encode()).hexdigest() + ".json"
path = Path(directory) / filename
require(path.is_file(), "Missing main result: " + filename)
record = json.loads(path.read_text())
require(record["key"] == item["key"], "Main result key mismatch")
for field in ("polynomial", "isolating_interval", "degree", "modified"):
require(record["argument"][field] == item[field], "Main result catalogue mismatch: " + field)
precisions, measurements = record["working_precisions"], record["measurements"]
require(len(precisions) == len(measurements) >= 2 and precisions == sorted(set(precisions)),
"Main result needs two ordered distinct precisions")
return record, path
def read_jsonl(path):
records = {}
for line_number, line in enumerate(Path(path).read_text().splitlines(), 1):
if not line.strip():
continue
record = json.loads(line)
kind = record.get("kind")
require(kind in ("ordinary", "modified", "quartic"), "Unknown JSONL kind on line " + str(line_number))
key = str(record["source_row"]) if kind == "quartic" else record["key"]
identity = kind + ":" + key
require(identity not in records, "Duplicate JSONL record: " + identity)
records[identity] = record
return records
def scalar_reading(item, kind, working, cache_dir=None):
fingerprint = object_hash({
"numerical_evaluator": NUMERICAL_FINGERPRINTS[kind],
"argument": item, "kind": kind, "working_digits": working,
"mpmath": mp.__version__, "python": platform.python_version(),
})
cache_path = Path(cache_dir) / (fingerprint + ".json") if cache_dir else None
if cache_path and cache_path.exists():
cached = json.loads(cache_path.read_text())
require(cached.get("fingerprint") == fingerprint, "Numerical cache fingerprint mismatch")
require(cached.get("reading_sha256") == object_hash(cached["reading"]), "Numerical cache content mismatch")
return {**cached["reading"], "cache_reused": True, "numerical_cache_fingerprint": fingerprint}
parameter, isolation = resolve_argument(item, working)
evaluate = cubic_jensen if kind == "ordinary" else modified_hypergeometric
result, evidence = evaluate(parameter, working)
reading = {"working_digits": working, "value": mp.nstr(result, working),
"isolation": isolation, **evidence}
if cache_path:
save_json(cache_path, {"fingerprint": fingerprint, "reading": reading,
"reading_sha256": object_hash(reading),
"numerical_evaluator": NUMERICAL_FINGERPRINTS[kind],
"loaded_evaluator_sha256": EVALUATOR_SHA256,
"status": "computed_not_yet_compared"})
return {**reading, "cache_reused": False, "numerical_cache_fingerprint": fingerprint}
def verify_scalar(item, kind, supplied, working_precisions, cache_dir=None):
readings = [scalar_reading(item, kind, working, cache_dir) for working in working_precisions]
with mp.workdps(max(working_precisions) + 40):
high = decimal(readings[-1]["value"])
low = decimal(readings[0]["value"])
zero = kind == "modified" and readings[-1]["isolation"]["exact_rational"] == "0"
comparisons = []
if zero:
require(high == low == 0, "Exact-zero computation failed")
require(supplied.get("exact") is True, "Modified zero requires explicit exact status")
require(all(decimal(value) == 0 for value in supplied["high_values"]), "Main modified zero is nonzero")
else:
comparisons.append(agreement(low, high, "Independent two-precision stability"))
for value in supplied["high_values"]:
comparisons.append(agreement(decimal(value), high, "Main versus independent " + kind))
if "stored" in supplied:
comparisons.append(check_stored(supplied["stored"], high, exact=zero))
controls = []
exact_parameter = readings[-1]["isolation"]["exact_rational"]
if kind == "ordinary" and exact_parameter in ("0", "27"):
target = dirichlet_control() * (3 if exact_parameter == "27" else 1)
controls.append({"control": "nonzero Dirichlet special value", **agreement(high, target, "N(0)/N(27) control")})
if exact_parameter in ("-216", "-27", "-8", "-1", "32", "54", "125") and kind == "ordinary" or exact_parameter == "8" and kind == "modified":
expected, evidence = elliptic_control(int(exact_parameter), kind, max(working_precisions))
controls.append({**evidence, **agreement(high, expected, "Proven elliptic L-value control")})
return {"status": "passed", "kind": kind, "key": item["key"],
"readings": readings, "comparisons": comparisons, "controls": controls}
def run_checks(options):
arguments, orbits = load_catalogue(options.catalogue)
records = read_jsonl(options.values) if options.values else {}
identities = ["ordinary:" + key for key in arguments]
identities += ["modified:" + key for key, item in arguments.items() if item["modified"]]
identities += ["quartic:" + row for row in orbits]
require(not records or set(records) <= set(identities), "JSONL contains unknown catalogue entries")
selected = options.select or identities
require(len(selected) == len(set(selected)) and set(selected) <= set(identities), "Invalid or repeated --select")
report = {"provenance": provenance(options.catalogue, options.values),
"working_precisions": options.working, "full_catalogue_requested": not options.select,
"selected": selected, "items": [], "status": "running"}
save_json(options.report, report)
output_dir = Path(options.report).resolve().parent / "per-key"
independent_values = {}
with deadline(options.timeout), mp.workdps(max(options.working) + 25):
control_value, control_evidence = cubic_jensen(mp.mpf(0), options.working[0])
report["nonzero_instrument_control"] = {
**control_evidence, **agreement(control_value, dirichlet_control(), "Pre-run nonzero Jensen control")}
save_json(options.report, report)
def verify_identity(identity):
kind, key = identity.split(":", 1)
dependency_results = []
if kind == "quartic":
orbit = orbits[key]
for argument_key in orbit["keys"]:
branch = "modified" if arguments[argument_key]["modified"] else "ordinary"
dependency_results.append(verify_identity(branch + ":" + argument_key))
require(identity in records, "Quartic output requires --values JSONL alongside --results")
source_record = records[identity]
input_hash = object_hash(source_record)
supplied = {"high_values": [source_record["value"]]}
elif options.results:
source_record, result_path = read_main_result(options.results, arguments[key])
input_hash = sha256(result_path)
supplied = {"high_values": [measurement[kind] for measurement in source_record["measurements"]],
"stored": source_record["values"][kind],
"exact": kind == "modified" and type(source_record["values"][kind]) is int and source_record["values"][kind] == 0}
if identity in records:
declared = records[identity]
supplied["high_values"].append(declared["value"])
if "stored" in declared:
supplied["stored"] = declared["stored"]
if kind == "modified" and supplied["exact"]:
require(declared.get("exact") is True, "Final JSONL must preserve exact-zero status")
input_hash = object_hash({"calculation_sha256": input_hash, "declared_jsonl_record": declared})
else:
require(identity in records, "Missing JSONL record " + identity)
source_record = records[identity]
input_hash = object_hash(source_record)
supplied = {"high_values": [source_record["value"]], "exact": source_record.get("exact", False)}
if "stored" in source_record:
supplied["stored"] = source_record["stored"]
fingerprint = object_hash({"evaluator": EVALUATOR_SHA256, "catalogue": sha256(options.catalogue),
"input": input_hash, "identity": identity, "working": options.working,
"source_hashes": report["provenance"]["hashes"],
"dependencies": [result["fingerprint"] for result in dependency_results]})
output_path = output_dir / (hashlib.sha256(identity.encode()).hexdigest() + ".json")
if output_path.exists() and not options.force:
existing = json.loads(output_path.read_text())
if existing.get("fingerprint") == fingerprint and existing.get("status") == "passed":
independent_values[identity] = existing
return existing
start = time.monotonic()
with deadline(options.timeout):
if kind != "quartic":
result = verify_scalar(arguments[key], kind, supplied, options.working,
cache_dir=HERE / "scalar-cache")
else:
with mp.workdps(max(options.working) + 25):
values = [decimal(result["readings"][-1]["value"]) for result in dependency_results]
determinant, evidence = signed_determinant(values, orbit["coefficients"], orbit["variant"])
require(mp.sign(determinant) == (1 if rational(orbit["r"]) > 0 else -1),
"Determinant sign disagrees with the source coefficient of L(E,2)")
low_values = [decimal(result["readings"][0]["value"]) for result in dependency_results]
low, low_evidence = signed_determinant(low_values, orbit["coefficients"], orbit["variant"])
comparisons = [agreement(determinant, low, "Determinant stability from independent low/high periods"),
agreement(determinant, decimal(source_record["value"]), "Main versus independent quartic")]
if "stored" in source_record:
comparisons.append(check_stored(source_record["stored"], determinant))
result = {"status": "passed", "kind": kind, "source_row": key,
"value": mp.nstr(determinant, max(options.working)), "evidence": evidence,
"low_expansion_check": low_evidence["expansion_check"], "comparisons": comparisons,
"r": str(orbit["r"]), "s": str(orbit["s"]),
"L_function_identity_numerically_checked": False}
result.update(fingerprint=fingerprint, seconds=round(time.monotonic() - start, 3), input_sha256=input_hash)
save_json(output_path, result)
independent_values[identity] = result
return result
for identity in selected:
report["current_identity"] = identity
report["current_started_utc"] = datetime.now(timezone.utc).isoformat()
save_json(options.report, report)
print(json.dumps({"identity": identity, "status": "started"}), flush=True)
try:
result = verify_identity(identity)
report["items"].append({"identity": identity, "status": result["status"], "fingerprint": result["fingerprint"]})
print(json.dumps({"identity": identity, "status": "passed"}), flush=True)
except Exception as error:
failure = {"identity": identity, "status": "failed", "error_type": type(error).__name__, "error": str(error)}
save_json(output_dir / (hashlib.sha256(identity.encode()).hexdigest() + ".failure.json"),
{**failure, "provenance": report["provenance"]})
report["items"].append(failure)
print(json.dumps(failure), flush=True)
save_json(options.report, report)
failures = [item for item in report["items"] if item["status"] != "passed"]
report.update(status="failed" if failures else "passed", full_catalogue_verified=not options.select and not failures,
failures=len(failures), passing=len(report["items"]) - len(failures), current_identity=None)
save_json(options.report, report)
return 1 if failures else 0
def run_smoke(options):
report = {"provenance": provenance(), "scope": "small documented controls, not full-catalogue verification",
"full_catalogue_verified": False, "items": [], "status": "running"}
save_json(options.report, report)
cases = [("ordinary", 0), ("ordinary", 27), ("ordinary", 32), ("modified", 0), ("modified", 8)]
for kind, parameter in cases:
start = time.monotonic()
try:
with deadline(options.timeout), mp.workdps(max(options.working) + 30):
if kind == "ordinary" and parameter in (0, 27):
target = dirichlet_control() * (3 if parameter == 27 else 1)
elif parameter == 0:
target = mp.mpf(0)
else:
target, _ = elliptic_control(parameter, kind, max(options.working) + 15)
item = {"key": str(parameter), "polynomial": [-parameter, 1], "degree": 1,
"isolating_interval": [str(parameter - 1), str(parameter + 1)],
"modified": -1 < parameter < 27}
supplied = {"high_values": [mp.nstr(target, max(options.working) + 15)], "exact": parameter == 0 and kind == "modified"}
result = verify_scalar(item, kind, supplied, options.working)
if parameter == 32:
independent_series = exterior_series(mp.mpf(32))
result["controls"].append(agreement(independent_series, target, "Exterior constant-term series at t=32"))
target_g4 = 3 * target / 4
measured_g4, evidence_g4 = g4_jensen(options.working[0])
result["controls"].append({"relation": "g(4)=3*N(32)/4=2*L'(E20,0)",
"value_g4": mp.nstr(measured_g4, options.working[0]),
"source": "Rogers-Zudilin 1012.3036v3 equation (2)",
**evidence_g4, **agreement(measured_g4, target_g4, "g(4) independent integral")})
result["seconds"] = round(time.monotonic() - start, 3)
report["items"].append(result)
print(json.dumps({"kind": kind, "key": str(parameter), "status": "passed", "seconds": result["seconds"]}), flush=True)
except Exception as error:
report["items"].append({"kind": kind, "key": str(parameter), "status": "failed",
"error_type": type(error).__name__, "error": str(error)})
print(json.dumps(report["items"][-1]), flush=True)
save_json(options.report, report)
report["status"] = "passed" if all(item["status"] == "passed" for item in report["items"]) else "failed"
save_json(options.report, report)
return 0 if report["status"] == "passed" else 1
def main():
parser = argparse.ArgumentParser(description=__doc__)
parser.add_argument("command", choices=("check", "smoke"))
parser.add_argument("--catalogue")
parser.add_argument("--results", help="Main per-key directory; filenames sha256(key).json")
parser.add_argument("--values", help="JSONL scalar records or supplementary quartic results")
parser.add_argument("--report", required=True)
parser.add_argument("--select", action="append", help="Explicit partial identity, e.g. ordinary:32")
parser.add_argument("--working", nargs=2, type=int, default=[125, 145])
parser.add_argument("--timeout", type=int, default=1800, help="Seconds per scalar/determinant item")
parser.add_argument("--force", action="store_true", help="Recompute even a matching passing checkpoint")
options = parser.parse_args()
require(125 <= options.working[0] < options.working[1], "Use distinct increasing working precisions >=125")
require(options.timeout > 0, "Timeout must be positive")
if options.command == "smoke":
return run_smoke(options)
require(options.catalogue and (options.values or options.results), "check requires catalogue and results/values")
try:
return run_checks(options)
except Exception as error:
save_json(options.report, {"status": "failed", "full_catalogue_verified": False,
"error_type": type(error).__name__, "error": str(error),
"provenance": provenance(options.catalogue, options.values)})
raise
NUMERICAL_FINGERPRINTS = {}
for numerical_kind, numerical_function in (("ordinary", cubic_jensen), ("modified", modified_hypergeometric)):
NUMERICAL_FINGERPRINTS[numerical_kind] = object_hash([
inspect.getsource(function) for function in (
rational, trim, polynomial_value, remainder, sturm_sequence,
variations, validate_argument, resolve_argument, numerical_function,
*([modified_boundary_series, dirichlet_control] if numerical_kind == "modified" else []),
)
])
if __name__ == "__main__":
try:
sys.exit(main())
except Exception as failure:
print(json.dumps({"status": "failed", "error_type": type(failure).__name__, "error": str(failure)}), file=sys.stderr)
sys.exit(2)