#!/usr/bin/env sage
"""
Exact algebraic hypotheses for the Problem 2.8 tail-lattice induction.

Run from the repository root with

    sage agent_outputs/tail_lattice/p28_lattice_hypotheses_certificate.sage

The script loads the independent all-N kernel certificate, then verifies:

  * J_N = diag(x,1,1,1) M_N is regular at x=0;
  * J_N(0) has the claimed rank-one factorization;
  * its image direction is the leading direction of H k_N;
  * the transformed 3F2 equation gives the exact row dependence.

The only non-machine step in the lattice closure is then the two-line DVR
lemma proved in TAIL_LATTICE_CLOSURE_REPORT.md.
"""

from sage.all import *
import os


HERE = os.path.dirname(os.path.abspath(__file__))
KERNEL = os.path.join(HERE, "p28_kernel_contiguity_certificate.sage")
if not os.path.exists(KERNEL):
    KERNEL = os.path.join(
        HERE, "..", "special_functions",
        "p28_kernel_contiguity_certificate.sage"
    )
load(KERNEL)


H = diagonal_matrix(K, [x, 1, 1, 1])
J = H*M


def value_at_zero(q):
    """Evaluate a simplified rational function at x=0."""
    q = K(q)
    numerator = q.numerator()
    denominator = q.denominator()
    value_denominator = denominator.subs({x: 0})
    assert value_denominator != 0
    return K(numerator.subs({x: 0}) / value_denominator)


J0 = Matrix(K, 4, 4, [value_at_zero(q) for q in J.list()])
a = 144*(u-1)**2 / (u*(3*u-2)*(3*u+2))
left = vector(K, [a, -1, -1, -1])
right = vector(K, [u**3, 3*u**2, 3*u, 1])

assert J0 == left.column()*right.row()
assert J0.rank() == 1
assert J0.column(0) != 0


# The x^(-1) coefficient of M controls the constant term of the transformed
# denominator.  It is another rank-one matrix, now with only its first row
# nonzero.
Mminus1 = Matrix(K, 4, 4, [
    value_at_zero(x*q) for q in M.list()
])
e0 = vector(K, [1, 0, 0, 0])
assert Mminus1 == e0.column()*(a*right).row()

# Therefore c_(N+1)/c_N is the first entry of a*right.
constant_ratio = a*u**3
assert constant_ratio == (
    144*(u-1)**2*u**2 / ((3*u-2)*(3*u+2))
)


# The first coefficient of
#
# 4F3(m,m+1/6,m+1/2,m+5/6;2m,2m,2m;z)
#
# is c.  Since z=-x+O(x^2), the leading direction of H*k_N is
# (1,-c,-c,-c)^T.
c = (
    m*(m+QQ(1)/6)*(m+QQ(1)/2)*(m+QQ(5)/6)
    / (2*m)**3
)
assert K(c - 1/a) == 0
tail_direction = vector(K, [1, -c, -c, -c])
assert left == a*tail_direction


# The row dependence is just the transformed 3F2 equation.  It is recorded
# here as a formal coefficient identity in the four symbols
# (y-1, theta*y, theta^2*y, theta^3*y).
Y = PolynomialRing(K, names=("f0", "f1", "f2", "f3"))
f0, f1, f2, f3 = Y.gens()
y = f0 + 1
ode = 72*f3 + 108*x*f2 + 46*x*f1 + 5*x*y

# C*k0 = -5/4 after B*k0=f.
Ck0 = 18*f3/x + QQ(5)/4*f0 + QQ(23)/2*f1 + 27*f2
assert Y(x*(Ck0 + QQ(5)/4) - ode/4) == 0

# Therefore 72 E3 + 108 x E2 + 46 x E1 + 5 x E0 = 0.
# The B-part cancels independently.
b0 = vector(K, [1, 0, 0, 0])
b1 = vector(K, [1, 1, 0, 0])
b2 = vector(K, [1, 2, 1, 0])
b3 = vector(K, [1, 3, 3, 1])
assert 72*b3 + 108*x*b2 + 46*x*b1 + 5*x*b0 == vector(
    K, [72 + 108*x + 46*x + 5*x,
        216 + 108*x + 46*x,
        216 + 108*x,
        72]
)

# In the full carrier the f*C contribution is killed by the ODE, and the
# displayed B combination equals -4*C after using the compact expression
# C=(18/x)b3+(5/4)b0+(23/2)b1+27b2.  Verify this exact cancellation.
Crow = 18*b3/x + QQ(5)/4*b0 + QQ(23)/2*b1 + 27*b2
Bcomb = 72*b3 + 108*x*b2 + 46*x*b1 + 5*x*b0
assert Bcomb == 4*x*Crow


print("PASS: exact rank-one/DVR hypotheses for all N")
print("J_N(0) = (a,-1,-1,-1)^T (u^3,3u^2,3u,1)")
print("a^{-1} is the first shifted 4F3 coefficient")