743 lines
19 KiB
Python
743 lines
19 KiB
Python
#!/usr/bin/env python3
|
|
"""Transparent exact equation checks for Ramanujan Challenge 2.8.
|
|
|
|
This verifier deliberately implements only:
|
|
|
|
* sparse multivariate polynomials over ``fractions.Fraction``;
|
|
* rational functions represented by numerator/denominator pairs;
|
|
* addition, multiplication, integer powers, and formal differentiation.
|
|
|
|
It does not call polynomial division, factorization, simplification, Gröbner
|
|
bases, a special-function library, a root finder, or numerical sampling.
|
|
Every obligation is entered as a displayed denominator-cleared identity and
|
|
passes only when every coefficient of the expanded numerator is exactly zero.
|
|
"""
|
|
|
|
from fractions import Fraction as F
|
|
|
|
|
|
VARIABLES = ("u", "x", "t", "j", "n", "k", "z")
|
|
NVARS = len(VARIABLES)
|
|
INDEX = {name: position for position, name in enumerate(VARIABLES)}
|
|
ZERO_EXPONENT = (0,) * NVARS
|
|
|
|
|
|
class Poly:
|
|
"""Sparse polynomial over QQ in the fixed variables above."""
|
|
|
|
def __init__(self, terms=None):
|
|
combined = {}
|
|
for exponent, coefficient in (terms or {}).items():
|
|
coefficient = F(coefficient)
|
|
if coefficient:
|
|
combined[tuple(exponent)] = (
|
|
combined.get(tuple(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("Poly 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 over the sparse polynomial ring."""
|
|
|
|
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))
|
|
|
|
|
|
def check_zero(label, expression):
|
|
expression = as_rat(expression)
|
|
assert expression.is_zero(), label
|
|
print("PASS:", label)
|
|
|
|
|
|
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))
|
|
]
|
|
|
|
|
|
def check_zero_matrix(label, matrix):
|
|
assert all(
|
|
as_rat(entry).is_zero()
|
|
for row in matrix
|
|
for entry in row
|
|
), label
|
|
print("PASS:", label)
|
|
|
|
|
|
u, x, t, j, n, k, z = [
|
|
Rat(Poly.variable(name)) for name in VARIABLES
|
|
]
|
|
SYMBOLS = dict(zip(VARIABLES, (u, x, t, j, n, k, z)))
|
|
|
|
|
|
def substitute_polynomial(polynomial, replacements):
|
|
"""Evaluate a sparse polynomial at rational-function 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):
|
|
"""Evaluate an unreduced rational function by cross multiplication."""
|
|
expression = as_rat(expression)
|
|
return (
|
|
substitute_polynomial(expression.numerator, replacements)
|
|
/ substitute_polynomial(expression.denominator, replacements)
|
|
)
|
|
|
|
|
|
def polynomial_coefficient(expression, name, degree):
|
|
"""Extract one 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, coefficient in expression.numerator.terms.items():
|
|
if exponent[position] == degree:
|
|
reduced = list(exponent)
|
|
reduced[position] = 0
|
|
terms[tuple(reduced)] = coefficient
|
|
return Rat(Poly(terms), expression.denominator)
|
|
|
|
|
|
# ---------------------------------------------------------------------------
|
|
# Four standalone tail/Ore factorizations.
|
|
# ---------------------------------------------------------------------------
|
|
|
|
m = (u - 1) / 2
|
|
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 + 58 + F(14, 9))
|
|
+ (-4*u**2 + 32*u - 63)
|
|
)
|
|
|
|
matrix = [
|
|
[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],
|
|
]
|
|
|
|
# The first column of the z-coordinate gauge
|
|
# -z M(2n+1,-z/(1-z)).
|
|
# It will be contracted below with four explicitly reconstructed horizontal
|
|
# row components to verify that the displayed d0+z*d1 step is induced by the
|
|
# authoritative matrix, rather than merely guessed and checked afterward.
|
|
matrix_nz = [
|
|
[
|
|
substitute_rational(
|
|
entry,
|
|
{"u": 2*n+1, "x": -z/(1-z)},
|
|
)
|
|
for entry in row
|
|
]
|
|
for row in matrix
|
|
]
|
|
gauge_first_column = [-z*matrix_nz[row][0] for row in range(4)]
|
|
|
|
|
|
def p_row(row, argument=t):
|
|
return sum(matrix[row][column] * argument**column for column in range(4))
|
|
|
|
|
|
P = [p_row(row) for row in range(4)]
|
|
|
|
|
|
def shifted_derivative(expression):
|
|
return (1-x)*x*expression.derivative("x") + (t+1)*expression
|
|
|
|
|
|
L_plus = (
|
|
(1-x)*t*(t+u)**3
|
|
+ x*(t+m+1)*(t+m+F(7, 6))*(t+m+F(3, 2))*(t+m+F(11, 6))
|
|
)
|
|
|
|
l0 = (u-1)*u*(3*u-2)*(3*u+2)*x/144
|
|
l1 = (
|
|
-576 + 864*u - 432*u**2 + 72*u**3
|
|
+ 580*x - 872*u*x + 405*u**2*x - 36*u**3*x
|
|
) / 72
|
|
l2 = (
|
|
432 - 432*u + 108*u**2
|
|
- 436*x + 405*u*x - 54*u**2*x
|
|
) / 36
|
|
l3 = (-12 + 6*u + 11*x - 2*u*x) / 2
|
|
|
|
D0 = shifted_derivative(P[0]) - P[1]
|
|
D1 = shifted_derivative(P[1]) - P[2]
|
|
D2 = shifted_derivative(P[2]) - P[3]
|
|
D3 = shifted_derivative(P[3]) + l3*P[3] + l2*P[2] + l1*P[1] + l0*P[0]
|
|
|
|
q3_numerator = (
|
|
-36 + 536*x - 297*u*x + 54*u**2*x
|
|
- 567*x**2 + 288*u*x**2 - 36*u**2*x**2
|
|
)
|
|
|
|
check_zero(
|
|
"cleared tail factorization D0=q0*L_plus",
|
|
u*(3*u-2)*(3*u+2)*x*D0 - 144*(u-1)**2*L_plus,
|
|
)
|
|
check_zero("cleared tail factorization D1=-L_plus", D1 + L_plus)
|
|
check_zero(
|
|
"cleared tail factorization D2=q2*L_plus",
|
|
2*D2 - (-2+7*x-2*u*x)*L_plus,
|
|
)
|
|
check_zero(
|
|
"cleared fourth companion closure D3=q3*L_plus",
|
|
36*D3 - q3_numerator*L_plus,
|
|
)
|
|
|
|
|
|
# Row-zero coefficient identities. These prove F_N=P_0(delta)F_(N+1)
|
|
# without any Ore division or finite sampling.
|
|
def A(argument):
|
|
return (
|
|
144*(u-1)**2*(argument+u)**3
|
|
/ (u*(3*u-2)*(3*u+2))
|
|
)
|
|
|
|
|
|
def P0(argument):
|
|
return p_row(0, argument)
|
|
|
|
|
|
def B(argument):
|
|
return P0(argument) - A(argument)/x
|
|
|
|
|
|
rho = -(3*u-2)*(3*u+2) / (144*(u-1)**2*u**2)
|
|
|
|
|
|
def coeff_cross_ratio(argument):
|
|
return (
|
|
(
|
|
(m+argument)
|
|
* (m+F(1, 6)+argument)
|
|
* (m+F(1, 2)+argument)
|
|
* (m+F(5, 6)+argument)
|
|
)
|
|
/ (m*(m+F(1, 6))*(m+F(1, 2))*(m+F(5, 6)))
|
|
* (
|
|
2*m*(2*m+1)
|
|
/ ((2*m+argument)*(2*m+argument+1))
|
|
)**3
|
|
)
|
|
|
|
|
|
def coeff_within_ratio(argument):
|
|
return (
|
|
(m+argument)
|
|
* (m+F(1, 6)+argument)
|
|
* (m+F(1, 2)+argument)
|
|
* (m+F(5, 6)+argument)
|
|
/ ((2*m+argument)**3*(argument+1))
|
|
)
|
|
|
|
|
|
check_zero("tail lowest-coefficient normalization", rho*(-A(0))-1)
|
|
check_zero(
|
|
"tail generic coefficient identity",
|
|
rho * (
|
|
-A(j)*coeff_cross_ratio(j)
|
|
+ (A(j-1)+B(j-1))
|
|
* coeff_cross_ratio(j-1)
|
|
/ coeff_within_ratio(j-1)
|
|
) - 1,
|
|
)
|
|
|
|
|
|
# ---------------------------------------------------------------------------
|
|
# Terminating denominator: explicit one-step coefficient induction.
|
|
# ---------------------------------------------------------------------------
|
|
|
|
normalization = 576*n**2*(2*n+1)**2 / ((6*n+1)*(6*n+5))
|
|
|
|
|
|
def d0(argument):
|
|
return (
|
|
72*(2*n+1)**2*(2*n-argument)**3
|
|
/ (n*(6*n+1)*(6*n+5))
|
|
)
|
|
|
|
|
|
def d1_polynomial(argument):
|
|
return (
|
|
5*argument - 51*argument**2 - 72*argument**3
|
|
- 5*n + 127*argument*n + 659*argument**2*n - 432*argument**3*n
|
|
- 76*n**2 - 1760*argument*n**2 + 3086*argument**2*n**2
|
|
- 864*argument**3*n**2
|
|
+ 1404*n**3 - 6536*argument*n**3 + 4500*argument**2*n**3
|
|
- 576*argument**3*n**3
|
|
+ 4360*n**4 - 7632*argument*n**4 + 2232*argument**2*n**4
|
|
+ 4320*n**5 - 3024*argument*n**5
|
|
+ 1440*n**6
|
|
)
|
|
|
|
|
|
def d1(argument):
|
|
return (
|
|
-d1_polynomial(argument)
|
|
/ (n*(2*n+1)*(6*n+1)*(6*n+5))
|
|
)
|
|
|
|
|
|
# Horizontal adjoint reconstruction. If p is the first row component, the
|
|
# other three are p_i(theta)p with the following explicit polynomials. Their
|
|
# contraction with the first column of -z M(2n+1,-z/(1-z)) must be exactly
|
|
# d0(theta)+z*d1(theta). This is the formerly implicit matrix-to-scalar
|
|
# bridge.
|
|
p_ops = [
|
|
Rat(1),
|
|
(
|
|
2*t*(
|
|
-36 + 216*n - 432*n**2
|
|
- 36*t + 216*n*t - 36*t**2
|
|
+ 23*z + 162*n*z + 216*n**2*z
|
|
- 54*t*z - 144*n*t*z + 36*t**2*z
|
|
)
|
|
/ (n*(1+2*n)*(1+6*n)*(5+6*n)*z)
|
|
),
|
|
(
|
|
36*t*(
|
|
4 - 12*n + 2*t + 3*z + 8*n*z - 2*t*z
|
|
)
|
|
/ (n*(1+2*n)*(1+6*n)*(5+6*n)*z)
|
|
),
|
|
(
|
|
72*t*(-1+z)
|
|
/ (n*(1+2*n)*(1+6*n)*(5+6*n)*z)
|
|
),
|
|
]
|
|
|
|
l_tail = (
|
|
t*(t+2*n-1)**3
|
|
- z*(t+n)*(t+n+F(1, 6))*(t+n+F(1, 2))*(t+n+F(5, 6))
|
|
)
|
|
l_coefficients = [
|
|
polynomial_coefficient(l_tail, "t", degree)
|
|
for degree in range(5)
|
|
]
|
|
|
|
|
|
def companion(parameter):
|
|
operator = (
|
|
t*(t+2*parameter-1)**3
|
|
- z*(t+parameter)
|
|
*(t+parameter+F(1, 6))
|
|
*(t+parameter+F(1, 2))
|
|
*(t+parameter+F(5, 6))
|
|
)
|
|
coefficients = [
|
|
polynomial_coefficient(operator, "t", degree)
|
|
for degree in range(5)
|
|
]
|
|
return [
|
|
[0, 1, 0, 0],
|
|
[0, 0, 1, 0],
|
|
[0, 0, 0, 1],
|
|
[
|
|
-coefficients[column]/coefficients[4]
|
|
for column in range(4)
|
|
],
|
|
]
|
|
|
|
|
|
def theta_operator(operator):
|
|
"""Left-coefficient Euler composition: theta Q = z Q_z + t Q."""
|
|
return z*operator.derivative("z") + t*operator
|
|
|
|
|
|
check_zero(
|
|
"horizontal reconstruction pi3",
|
|
p_ops[3] - l_coefficients[4]/l_coefficients[0]*t,
|
|
)
|
|
check_zero(
|
|
"horizontal reconstruction pi2",
|
|
p_ops[2]
|
|
- l_coefficients[3]/l_coefficients[4]*p_ops[3]
|
|
+ theta_operator(p_ops[3]),
|
|
)
|
|
check_zero(
|
|
"horizontal reconstruction pi1",
|
|
p_ops[1]
|
|
- l_coefficients[2]/l_coefficients[4]*p_ops[3]
|
|
+ theta_operator(p_ops[2]),
|
|
)
|
|
|
|
terminating_operator = (
|
|
t*(t-2*n)**3
|
|
- z*(t-n)*(t-n-F(1, 6))*(t-n-F(1, 2))*(t-n-F(5, 6))
|
|
)
|
|
adjoint_factor = -72/(n*(2*n+1)*(6*n+1)*(6*n+5)*z)
|
|
check_zero(
|
|
"horizontal reconstruction closes to the terminating operator",
|
|
theta_operator(p_ops[1]) + 1
|
|
- l_coefficients[1]/l_coefficients[4]*p_ops[3]
|
|
- adjoint_factor*terminating_operator,
|
|
)
|
|
|
|
gauge = [[-z*entry for entry in row] for row in matrix_nz]
|
|
left_gauge = matrix_multiply(companion(n), gauge)
|
|
right_gauge = matrix_multiply(gauge, companion(n+1))
|
|
gauge_residual = [
|
|
[
|
|
left_gauge[row][column]
|
|
- z*gauge[row][column].derivative("z")
|
|
- right_gauge[row][column]
|
|
for column in range(4)
|
|
]
|
|
for row in range(4)
|
|
]
|
|
check_zero_matrix(
|
|
"all sixteen authoritative differential-gauge equations",
|
|
gauge_residual,
|
|
)
|
|
|
|
matrix_induced_step = sum(
|
|
p_ops[row] * gauge_first_column[row]
|
|
for row in range(4)
|
|
)
|
|
check_zero(
|
|
"authoritative matrix induces the displayed d0+z*d1 scalar step",
|
|
matrix_induced_step - d0(t) - z*d1(t),
|
|
)
|
|
|
|
# Base polynomial from the compact denominator row:
|
|
# q_0(x)=18/x+159/4, Q_0=x q_0, p_1=(1-z)Q_0(-z/(1-z)).
|
|
q_zero = 18/x + F(159, 4)
|
|
Q_zero = x*q_zero
|
|
p_one = (1-z)*substitute_rational(Q_zero, {"x": -z/(1-z)})
|
|
check_zero(
|
|
"terminating base polynomial from the compact denominator row",
|
|
p_one - 18*(1-F(77, 24)*z),
|
|
)
|
|
|
|
|
|
def apply_euler_operator(operator, function, maximum_degree=3):
|
|
result = Rat(0)
|
|
theta_power = function
|
|
for degree in range(maximum_degree+1):
|
|
result += polynomial_coefficient(operator, "t", degree)*theta_power
|
|
theta_power = z*theta_power.derivative("z")
|
|
return result
|
|
|
|
|
|
compact_denominator = [
|
|
18/x + F(159, 4),
|
|
54/x + F(131, 2),
|
|
54/x + 27,
|
|
18/x,
|
|
]
|
|
base_horizontal_row = [
|
|
-z*substitute_rational(entry, {"x": -z/(1-z)})
|
|
for entry in compact_denominator
|
|
]
|
|
for row in range(4):
|
|
reconstructed = apply_euler_operator(
|
|
substitute_rational(p_ops[row], {"n": 1}),
|
|
p_one,
|
|
)
|
|
check_zero(
|
|
f"base horizontal-row reconstruction component {row}",
|
|
base_horizontal_row[row] - reconstructed,
|
|
)
|
|
|
|
base_terminating_operator = substitute_rational(
|
|
terminating_operator,
|
|
{"n": 1},
|
|
)
|
|
check_zero(
|
|
"base polynomial satisfies the terminating operator",
|
|
apply_euler_operator(
|
|
base_terminating_operator,
|
|
p_one,
|
|
maximum_degree=4,
|
|
),
|
|
)
|
|
|
|
base_times_companion = matrix_multiply(
|
|
[base_horizontal_row],
|
|
companion(Rat(1)),
|
|
)[0]
|
|
for column in range(4):
|
|
check_zero(
|
|
f"base horizontal adjoint residual component {column}",
|
|
z*base_horizontal_row[column].derivative("z")
|
|
+ base_times_companion[column],
|
|
)
|
|
|
|
|
|
within_ratio = (
|
|
(k-1-n)
|
|
* (k-1-n-F(1, 6))
|
|
* (k-1-n-F(1, 2))
|
|
* (k-1-n-F(5, 6))
|
|
/ ((k-2*n)**3*k)
|
|
)
|
|
cross_ratio = (
|
|
(n+1)/(n+1-k)
|
|
* (n+1+F(1, 6))/(n+1+F(1, 6)-k)
|
|
* (n+1+F(1, 2))/(n+1+F(1, 2)-k)
|
|
* (n+1+F(5, 6))/(n+1+F(5, 6)-k)
|
|
* ((k-2*n-1)*(k-2*n)/(2*n*(2*n+1)))**3
|
|
)
|
|
top_ratio = (
|
|
-(n+F(7, 6))*(n+F(3, 2))*(n+F(11, 6))
|
|
/ (8*(2*n+1)**3)
|
|
)
|
|
|
|
check_zero("terminating constant-term normalization", d0(0)-normalization)
|
|
check_zero(
|
|
"terminating generic coefficient induction",
|
|
d0(k) + d1(k-1)/within_ratio - normalization*cross_ratio,
|
|
)
|
|
check_zero(
|
|
"terminating top-degree boundary",
|
|
d1(n) - normalization*top_ratio,
|
|
)
|
|
|
|
|
|
# ---------------------------------------------------------------------------
|
|
# Two formerly named special-function steps, reduced to coefficients.
|
|
# ---------------------------------------------------------------------------
|
|
|
|
ascension_left_ratio = (
|
|
(k-F(5, 6))*(k-F(1, 2))*(k-F(1, 6))/k**3
|
|
)
|
|
ascension_right_ratio = (
|
|
(k-1)
|
|
* (k-F(5, 6))
|
|
* (k-F(1, 2))
|
|
* (k-F(1, 6))
|
|
/ (k**3*(k-1))
|
|
)
|
|
check_zero("ascension base coefficient", F(1, 6)*F(1, 2)*F(5, 6)-F(5, 72))
|
|
check_zero(
|
|
"ascension consecutive-coefficient ratio",
|
|
ascension_left_ratio-ascension_right_ratio,
|
|
)
|
|
|
|
theta_product = (
|
|
(t+F(1, 6))*(t+F(1, 2))*(t+F(5, 6))
|
|
)
|
|
theta_product_expanded = (
|
|
t**3 + F(3, 2)*t**2 + F(23, 36)*t + F(5, 72)
|
|
)
|
|
check_zero(
|
|
"3F2 Euler-operator coefficient expansion",
|
|
theta_product-theta_product_expanded,
|
|
)
|
|
|
|
|
|
print("PASS: all standalone exact-equation obligations")
|
|
print("No division algorithm, CAS simplifier, root finder, or sampling was used.")
|