allaun
1d5273264d
- Positive-cone transport certificates (p28_positive_cone.py, POSITIVE_CONE_CERTIFICATE.md, POSITIVE_CONE_MANUSCRIPT_SECTION.tex) - Optimized differential gauge (p28_optimized_gauge.py, OPTIMIZED_GAUGE_CERTIFICATE.md) - Adversarial provenance supplements (p28_mutation_sensitivity.py, solution_pre_positive_cone.tex) - FAMM SCARS advisory records (FAMM_SCARS.md, p28_famm_scars.json, p28_famm_scars_validator.py) - Overview documentation (OPTIONAL_IMPROVEMENTS.md, ADVERSARIAL_AUDIT.md) These are independent, replayable supplements developed after the original exact closure. They can be verified independently with 'bash run_checks.sh' in the certificates directory.
634 lines
17 KiB
Python
634 lines
17 KiB
Python
#!/usr/bin/env python3
|
|
"""Optimized exact gauge certificate for Ramanujan Challenge Problem 2.8.
|
|
|
|
The existing all-purpose standalone checker intentionally leaves rational
|
|
functions unreduced. That is maximally transparent, but the direct
|
|
sixteen-entry differential-gauge calculation creates very large temporary
|
|
denominators.
|
|
|
|
This independent checker first proves the exact decomposition
|
|
|
|
D M(u,x) = J0(u) + x J1(u) + x^2 J2(u),
|
|
D = diag(x,1,1,1),
|
|
|
|
and the additional relation
|
|
|
|
J2 = e4 * ((2u-9)/2) * b.
|
|
|
|
It then clears the common z-denominators before forming the gauge residual.
|
|
Every assertion is an equality in a sparse polynomial ring over QQ. The
|
|
implementation provides only addition, multiplication, integer powers,
|
|
formal differentiation, substitution, and coefficient extraction. It does
|
|
not call a simplifier, polynomial division, factorizer, Groebner basis,
|
|
special-function library, root finder, or numerical sampler.
|
|
"""
|
|
|
|
from fractions import Fraction as F
|
|
from time import perf_counter
|
|
|
|
|
|
START_TIME = perf_counter()
|
|
|
|
VARIABLES = ("u", "x", "n", "z", "t")
|
|
NVARS = len(VARIABLES)
|
|
INDEX = {name: position for position, name in enumerate(VARIABLES)}
|
|
ZERO_EXPONENT = (0,) * NVARS
|
|
|
|
|
|
class Poly:
|
|
"""Sparse multivariate polynomial over QQ."""
|
|
|
|
def __init__(self, terms=None):
|
|
combined = {}
|
|
for exponent, coefficient in (terms or {}).items():
|
|
exponent = tuple(exponent)
|
|
coefficient = F(coefficient)
|
|
if coefficient:
|
|
combined[exponent] = (
|
|
combined.get(exponent, F(0)) + coefficient
|
|
)
|
|
self.terms = {
|
|
exponent: coefficient
|
|
for exponent, coefficient in combined.items()
|
|
if coefficient
|
|
}
|
|
|
|
@staticmethod
|
|
def constant(value):
|
|
value = F(value)
|
|
return Poly({ZERO_EXPONENT: value}) if value else Poly()
|
|
|
|
@staticmethod
|
|
def variable(name):
|
|
exponent = [0] * NVARS
|
|
exponent[INDEX[name]] = 1
|
|
return Poly({tuple(exponent): F(1)})
|
|
|
|
def __add__(self, other):
|
|
other = as_poly(other)
|
|
terms = dict(self.terms)
|
|
for exponent, coefficient in other.terms.items():
|
|
terms[exponent] = (
|
|
terms.get(exponent, F(0)) + coefficient
|
|
)
|
|
return Poly(terms)
|
|
|
|
__radd__ = __add__
|
|
|
|
def __neg__(self):
|
|
return Poly({
|
|
exponent: -coefficient
|
|
for exponent, coefficient in self.terms.items()
|
|
})
|
|
|
|
def __sub__(self, other):
|
|
return self + (-as_poly(other))
|
|
|
|
def __rsub__(self, other):
|
|
return as_poly(other) - self
|
|
|
|
def __mul__(self, other):
|
|
other = as_poly(other)
|
|
terms = {}
|
|
for left_exp, left_coefficient in self.terms.items():
|
|
for right_exp, right_coefficient in other.terms.items():
|
|
exponent = tuple(
|
|
left_exp[position] + right_exp[position]
|
|
for position in range(NVARS)
|
|
)
|
|
terms[exponent] = (
|
|
terms.get(exponent, F(0))
|
|
+ left_coefficient * right_coefficient
|
|
)
|
|
return Poly(terms)
|
|
|
|
__rmul__ = __mul__
|
|
|
|
def __pow__(self, exponent):
|
|
if exponent < 0:
|
|
raise ValueError("polynomial powers must be nonnegative")
|
|
result = Poly.constant(1)
|
|
base = self
|
|
power = exponent
|
|
while power:
|
|
if power & 1:
|
|
result = result * base
|
|
base = base * base
|
|
power //= 2
|
|
return result
|
|
|
|
def derivative(self, name):
|
|
position = INDEX[name]
|
|
terms = {}
|
|
for exponent, coefficient in self.terms.items():
|
|
degree = exponent[position]
|
|
if degree:
|
|
new_exponent = list(exponent)
|
|
new_exponent[position] -= 1
|
|
terms[tuple(new_exponent)] = coefficient * degree
|
|
return Poly(terms)
|
|
|
|
def is_zero(self):
|
|
return not self.terms
|
|
|
|
|
|
def as_poly(value):
|
|
if isinstance(value, Poly):
|
|
return value
|
|
return Poly.constant(value)
|
|
|
|
|
|
class Rat:
|
|
"""Unreduced rational function represented by two sparse polynomials."""
|
|
|
|
def __init__(self, numerator=0, denominator=1):
|
|
self.numerator = as_poly(numerator)
|
|
self.denominator = as_poly(denominator)
|
|
if self.denominator.is_zero():
|
|
raise ZeroDivisionError("zero polynomial denominator")
|
|
|
|
def __add__(self, other):
|
|
other = as_rat(other)
|
|
return Rat(
|
|
self.numerator * other.denominator
|
|
+ other.numerator * self.denominator,
|
|
self.denominator * other.denominator,
|
|
)
|
|
|
|
__radd__ = __add__
|
|
|
|
def __neg__(self):
|
|
return Rat(-self.numerator, self.denominator)
|
|
|
|
def __sub__(self, other):
|
|
return self + (-as_rat(other))
|
|
|
|
def __rsub__(self, other):
|
|
return as_rat(other) - self
|
|
|
|
def __mul__(self, other):
|
|
other = as_rat(other)
|
|
return Rat(
|
|
self.numerator * other.numerator,
|
|
self.denominator * other.denominator,
|
|
)
|
|
|
|
__rmul__ = __mul__
|
|
|
|
def __truediv__(self, other):
|
|
other = as_rat(other)
|
|
if other.numerator.is_zero():
|
|
raise ZeroDivisionError("division by the zero rational function")
|
|
return Rat(
|
|
self.numerator * other.denominator,
|
|
self.denominator * other.numerator,
|
|
)
|
|
|
|
def __rtruediv__(self, other):
|
|
return as_rat(other) / self
|
|
|
|
def __pow__(self, exponent):
|
|
if exponent >= 0:
|
|
return Rat(
|
|
self.numerator ** exponent,
|
|
self.denominator ** exponent,
|
|
)
|
|
return Rat(
|
|
self.denominator ** (-exponent),
|
|
self.numerator ** (-exponent),
|
|
)
|
|
|
|
def derivative(self, name):
|
|
return Rat(
|
|
self.numerator.derivative(name) * self.denominator
|
|
- self.numerator * self.denominator.derivative(name),
|
|
self.denominator ** 2,
|
|
)
|
|
|
|
def is_zero(self):
|
|
return self.numerator.is_zero()
|
|
|
|
|
|
def as_rat(value):
|
|
if isinstance(value, Rat):
|
|
return value
|
|
if isinstance(value, Poly):
|
|
return Rat(value)
|
|
return Rat(F(value))
|
|
|
|
|
|
u, x, n, z, t = [
|
|
Rat(Poly.variable(name)) for name in VARIABLES
|
|
]
|
|
SYMBOLS = dict(zip(VARIABLES, (u, x, n, z, t)))
|
|
|
|
|
|
def substitute_polynomial(polynomial, replacements):
|
|
result = Rat(0)
|
|
for exponent, coefficient in polynomial.terms.items():
|
|
term = Rat(coefficient)
|
|
for position, degree in enumerate(exponent):
|
|
if degree:
|
|
name = VARIABLES[position]
|
|
term *= replacements.get(name, SYMBOLS[name]) ** degree
|
|
result += term
|
|
return result
|
|
|
|
|
|
def substitute_rational(expression, replacements):
|
|
expression = as_rat(expression)
|
|
return (
|
|
substitute_polynomial(expression.numerator, replacements)
|
|
/ substitute_polynomial(expression.denominator, replacements)
|
|
)
|
|
|
|
|
|
def coefficient(expression, name, degree):
|
|
"""Extract a coefficient when the denominator omits ``name``."""
|
|
expression = as_rat(expression)
|
|
position = INDEX[name]
|
|
assert all(
|
|
exponent[position] == 0
|
|
for exponent in expression.denominator.terms
|
|
)
|
|
terms = {}
|
|
for exponent, value in expression.numerator.terms.items():
|
|
if exponent[position] == degree:
|
|
reduced = list(exponent)
|
|
reduced[position] = 0
|
|
terms[tuple(reduced)] = value
|
|
return Rat(Poly(terms), expression.denominator)
|
|
|
|
|
|
def matrix_multiply(left, right):
|
|
return [
|
|
[
|
|
sum(
|
|
left[row][middle] * right[middle][column]
|
|
for middle in range(len(right))
|
|
)
|
|
for column in range(len(right[0]))
|
|
]
|
|
for row in range(len(left))
|
|
]
|
|
|
|
|
|
OBLIGATIONS = 0
|
|
|
|
|
|
def check_zero(label, expression):
|
|
global OBLIGATIONS
|
|
assert as_rat(expression).is_zero(), label
|
|
OBLIGATIONS += 1
|
|
|
|
|
|
def check_matrix_entries(label, matrix):
|
|
for row in range(len(matrix)):
|
|
for column in range(len(matrix[0])):
|
|
check_zero(
|
|
f"{label}, entry ({row + 1},{column + 1})",
|
|
matrix[row][column],
|
|
)
|
|
print(f"PASS: {label} ({len(matrix) * len(matrix[0])} entries)")
|
|
|
|
|
|
def matrix_subtract(left, right):
|
|
return [
|
|
[
|
|
left[row][column] - right[row][column]
|
|
for column in range(len(left[0]))
|
|
]
|
|
for row in range(len(left))
|
|
]
|
|
|
|
|
|
# ---------------------------------------------------------------------------
|
|
# Authoritative matrix and the exact J0+xJ1+x^2J2 decomposition.
|
|
# ---------------------------------------------------------------------------
|
|
|
|
r = 1 / x
|
|
w = u * (3*u - 2) * (3*u + 2)
|
|
|
|
a1 = (
|
|
r * (144*u**5 - 288*u**4 + 144*u**3)
|
|
- 99*u**5 + 333*u**4 - 229*u**3 - 114*u**2 + 40*u + 64
|
|
)
|
|
a2 = (
|
|
r * (432*u**4 - 864*u**3 + 432*u**2)
|
|
- 243*u**4 + 909*u**3 - 868*u**2 - 80*u + 272
|
|
)
|
|
a3 = (
|
|
r * (432*u**3 - 864*u**2 + 432*u)
|
|
- 153*u**3 + 648*u**2 - 860*u + 360
|
|
)
|
|
a4 = r * 144 * (u - 1)**2
|
|
|
|
b1 = (
|
|
r * (-144*u**3)
|
|
+ 9*u**4 + 63*u**3 + 158*u**2 + 168*u + 64
|
|
)
|
|
b2 = (
|
|
r * (216*u**2)
|
|
+ 36*u**3 - 189*u**2 - 316*u - 168
|
|
)
|
|
b3 = (
|
|
r * (108*u)
|
|
+ 54*u**2 - 189*u - 158
|
|
)
|
|
|
|
c1 = (
|
|
r**2 * (-288*u**3)
|
|
+ r * (54*u**4 + 378*u**3 + 948*u**2 + 1008*u + 384)
|
|
+ 18*u**5 + 45*u**4 - 251*u**3 - 1086*u**2 - 1384*u - 576
|
|
)
|
|
c2 = (
|
|
r**2 * (-432*u**2)
|
|
+ r * (153*u**4 - 657*u**3 + 1292*u**2 + 2064*u + 1072)
|
|
- 72*u**4 + 702*u**3 - 1069*u**2 - 2508*u - 1512
|
|
)
|
|
c3 = (
|
|
r**2 * (-216*u)
|
|
+ r * (180*u**3 - 891*u**2 + 1450*u + 1116)
|
|
- 108*u**3 + 864*u**2 - 1385*u - 1422
|
|
)
|
|
c4 = (
|
|
r**2 * (-4)
|
|
+ r * (6*u**2 - 33*u + F(536, 9))
|
|
- 4*u**2 + 32*u - 63
|
|
)
|
|
|
|
matrix_m = [
|
|
[a1/w, a2/w, a3/w, a4/w],
|
|
[-u**3, -3*u**2, -3*u, -1],
|
|
[
|
|
x*b1/144,
|
|
-x*b2/72,
|
|
-x*b3/36,
|
|
x*(-2*r - (2*u - 7))/2,
|
|
],
|
|
[x**2*c1/288, x**2*c2/144, x**2*c3/72, x**2*c4/4],
|
|
]
|
|
|
|
v = [u**3, 3*u**2, 3*u, 1]
|
|
alpha = 144*(u - 1)**2 / w
|
|
j0 = (
|
|
[[alpha*v[column] for column in range(4)]]
|
|
+ [[-v[column] for column in range(4)] for _ in range(3)]
|
|
)
|
|
|
|
a_finite = [
|
|
-99*u**5 + 333*u**4 - 229*u**3 - 114*u**2 + 40*u + 64,
|
|
-243*u**4 + 909*u**3 - 868*u**2 - 80*u + 272,
|
|
-153*u**3 + 648*u**2 - 860*u + 360,
|
|
0,
|
|
]
|
|
b_row = [
|
|
(u + 1)*(u + 2)*(3*u + 4)*(3*u + 8)/144,
|
|
(-36*u**3 + 189*u**2 + 316*u + 168)/72,
|
|
(-54*u**2 + 189*u + 158)/36,
|
|
(7 - 2*u)/2,
|
|
]
|
|
c_row = [
|
|
(u + 1)*(u + 2)*(3*u + 4)*(3*u + 8)/48,
|
|
(153*u**4 - 657*u**3 + 1292*u**2 + 2064*u + 1072)/144,
|
|
(180*u**3 - 891*u**2 + 1450*u + 1116)/72,
|
|
(54*u**2 - 297*u + 536)/36,
|
|
]
|
|
j1 = [
|
|
[entry/w for entry in a_finite],
|
|
[0, 0, 0, 0],
|
|
b_row,
|
|
c_row,
|
|
]
|
|
|
|
# J2 is entered independently from the finite c_i terms. The subsequent
|
|
# check against e4*beta*b is therefore not true by construction.
|
|
j2 = [
|
|
[0, 0, 0, 0],
|
|
[0, 0, 0, 0],
|
|
[0, 0, 0, 0],
|
|
[
|
|
(
|
|
18*u**5 + 45*u**4 - 251*u**3
|
|
- 1086*u**2 - 1384*u - 576
|
|
)/288,
|
|
(
|
|
-72*u**4 + 702*u**3 - 1069*u**2
|
|
- 2508*u - 1512
|
|
)/144,
|
|
(-108*u**3 + 864*u**2 - 1385*u - 1422)/72,
|
|
(-4*u**2 + 32*u - 63)/4,
|
|
],
|
|
]
|
|
|
|
beta = (2*u - 9)/2
|
|
j2_rank_one = [
|
|
[0, 0, 0, 0],
|
|
[0, 0, 0, 0],
|
|
[0, 0, 0, 0],
|
|
[beta*entry for entry in b_row],
|
|
]
|
|
check_matrix_entries(
|
|
"J2=e4*beta*b proportionality",
|
|
matrix_subtract(j2, j2_rank_one),
|
|
)
|
|
|
|
d_times_m = [
|
|
[
|
|
x*matrix_m[row][column] if row == 0
|
|
else matrix_m[row][column]
|
|
for column in range(4)
|
|
]
|
|
for row in range(4)
|
|
]
|
|
j_decomposition = [
|
|
[
|
|
j0[row][column]
|
|
+ x*j1[row][column]
|
|
+ x**2*j2[row][column]
|
|
for column in range(4)
|
|
]
|
|
for row in range(4)
|
|
]
|
|
check_matrix_entries(
|
|
"D*M=J0+x*J1+x^2*J2",
|
|
matrix_subtract(d_times_m, j_decomposition),
|
|
)
|
|
|
|
|
|
# ---------------------------------------------------------------------------
|
|
# Common-denominator clearing after x=-z/(1-z), u=2n+1.
|
|
# ---------------------------------------------------------------------------
|
|
|
|
substitutions = {"u": 2*n + 1}
|
|
j0_n = [
|
|
[substitute_rational(entry, substitutions) for entry in row]
|
|
for row in j0
|
|
]
|
|
j1_n = [
|
|
[substitute_rational(entry, substitutions) for entry in row]
|
|
for row in j1
|
|
]
|
|
j2_n = [
|
|
[substitute_rational(entry, substitutions) for entry in row]
|
|
for row in j2
|
|
]
|
|
|
|
e_diagonal = [1 - z, -z, -z, -z]
|
|
d_z = (1 - z)**2
|
|
g_bar = [
|
|
[
|
|
e_diagonal[row] * (
|
|
(1 - z)**2*j0_n[row][column]
|
|
- z*(1 - z)*j1_n[row][column]
|
|
+ z**2*j2_n[row][column]
|
|
)
|
|
for column in range(4)
|
|
]
|
|
for row in range(4)
|
|
]
|
|
|
|
# Bind the cleared formula directly to the authoritative matrix, rather than
|
|
# relying only on the already-checked decomposition.
|
|
matrix_nz = [
|
|
[
|
|
substitute_rational(
|
|
entry,
|
|
{"u": 2*n + 1, "x": -z/(1 - z)},
|
|
)
|
|
for entry in row
|
|
]
|
|
for row in matrix_m
|
|
]
|
|
g_direct = [[-z*entry for entry in row] for row in matrix_nz]
|
|
check_matrix_entries(
|
|
"Gbar=(1-z)^2*(-z*M) after the exact substitution",
|
|
[
|
|
[
|
|
g_bar[row][column] - d_z*g_direct[row][column]
|
|
for column in range(4)
|
|
]
|
|
for row in range(4)
|
|
],
|
|
)
|
|
|
|
|
|
def tail_operator(parameter):
|
|
return (
|
|
t*(t + 2*parameter - 1)**3
|
|
- z*(t + parameter)
|
|
*(t + parameter + F(1, 6))
|
|
*(t + parameter + F(1, 2))
|
|
*(t + parameter + F(5, 6))
|
|
)
|
|
|
|
|
|
def companion_and_cleared(parameter):
|
|
coefficients = [
|
|
coefficient(tail_operator(parameter), "t", degree)
|
|
for degree in range(5)
|
|
]
|
|
companion = [
|
|
[0, 1, 0, 0],
|
|
[0, 0, 1, 0],
|
|
[0, 0, 0, 1],
|
|
[-coefficients[column]/coefficients[4] for column in range(4)],
|
|
]
|
|
cleared = [
|
|
[0, 1 - z, 0, 0],
|
|
[0, 0, 1 - z, 0],
|
|
[0, 0, 0, 1 - z],
|
|
[-coefficients[column] for column in range(4)],
|
|
]
|
|
return companion, cleared
|
|
|
|
|
|
companion_n, c_bar_n = companion_and_cleared(n)
|
|
companion_n1, c_bar_n1 = companion_and_cleared(n + 1)
|
|
|
|
check_matrix_entries(
|
|
"Cbar_n=(1-z)*C_n",
|
|
[
|
|
[
|
|
c_bar_n[row][column]
|
|
- (1 - z)*companion_n[row][column]
|
|
for column in range(4)
|
|
]
|
|
for row in range(4)
|
|
],
|
|
)
|
|
check_matrix_entries(
|
|
"Cbar_(n+1)=(1-z)*C_(n+1)",
|
|
[
|
|
[
|
|
c_bar_n1[row][column]
|
|
- (1 - z)*companion_n1[row][column]
|
|
for column in range(4)
|
|
]
|
|
for row in range(4)
|
|
],
|
|
)
|
|
|
|
# Verify the quotient-rule conversion on every actual Gbar entry:
|
|
#
|
|
# d(1-z) theta(Gbar/d)
|
|
# = (1-z) z Gbar' + 2z Gbar, d=(1-z)^2.
|
|
quotient_rule_residual = []
|
|
for row in range(4):
|
|
residual_row = []
|
|
for column in range(4):
|
|
rational_entry = g_bar[row][column] / d_z
|
|
left = d_z*(1 - z)*z*rational_entry.derivative("z")
|
|
right = (
|
|
(1 - z)*z*g_bar[row][column].derivative("z")
|
|
+ 2*z*g_bar[row][column]
|
|
)
|
|
residual_row.append(left - right)
|
|
quotient_rule_residual.append(residual_row)
|
|
check_matrix_entries(
|
|
"entrywise quotient-rule clearing",
|
|
quotient_rule_residual,
|
|
)
|
|
|
|
|
|
# ---------------------------------------------------------------------------
|
|
# Sixteen-entry cleared polynomial gauge, checked coefficient by coefficient.
|
|
# ---------------------------------------------------------------------------
|
|
|
|
left_gauge = matrix_multiply(c_bar_n, g_bar)
|
|
right_gauge = matrix_multiply(g_bar, c_bar_n1)
|
|
cleared_residual = [
|
|
[
|
|
left_gauge[row][column]
|
|
- (1 - z)*z*g_bar[row][column].derivative("z")
|
|
- 2*z*g_bar[row][column]
|
|
- right_gauge[row][column]
|
|
for column in range(4)
|
|
]
|
|
for row in range(4)
|
|
]
|
|
|
|
# Gbar has z-degree at most three and Cbar has z-degree at most one.
|
|
# Therefore every cleared residual has z-degree at most four. Checking all
|
|
# five coefficients of all sixteen entries is a complete polynomial check.
|
|
for row in range(4):
|
|
for column in range(4):
|
|
for degree in range(5):
|
|
check_zero(
|
|
(
|
|
"cleared gauge coefficient "
|
|
f"entry ({row + 1},{column + 1}), z^{degree}"
|
|
),
|
|
coefficient(cleared_residual[row][column], "z", degree),
|
|
)
|
|
print(
|
|
"PASS: cleared gauge entry "
|
|
f"({row + 1},{column + 1}) coefficients z^0,...,z^4"
|
|
)
|
|
|
|
|
|
ELAPSED = perf_counter() - START_TIME
|
|
print("PASS: optimized denominator-cleared differential gauge")
|
|
print(f"PASS: {OBLIGATIONS} exact scalar obligations")
|
|
print(f"Runtime: {ELAPSED:.6f} seconds")
|
|
print("No simplifier, division algorithm, factorizer, root finder, or sampling.")
|