#!/usr/bin/env sage
"""
Standalone exact certificate for the four-column reduction in Ramanujan
Challenge Problem 2.8.

It certifies the algebraic part of the cyclic-frame lemma:

* the exact balanced limit S;
* charpoly(S)=Q_R/R^2;
* an explicit left eigenvector w_rho;
* every coordinate of w_rho is nonzero; and
* the limiting cyclic frame [e1,S e1,S^2 e1,S^3 e1] is invertible.

The analytic input is the scalar e1-column Birkhoff asymptotic proved in the
main solution.  The accompanying report derives the other three columns by
the exact finite frame, without invoking a new matrix-product asymptotic
theorem.
"""
from sage.all import *

Pn.<n> = PolynomialRing(QQ)
Fn = Pn.fraction_field()
R = QQ(151931373056001)


def authoritative_matrix(u, R):
    """Problem 2.8 transfer after 236337691420383=(14R-567)/9."""
    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+QQ(14)/9) \
        + (-4*u^2+32*u-63)

    return matrix(Fn, [
        [a1/w, a2/w, a3/w, a4/w],
        [-u^3, -3*u^2, -3*u, -1],
        [b1/(144*R), -b2/(72*R), -b3/(36*R),
         (-2*R-(2*u-7))/(2*R)],
        [c1/(288*R^2), c2/(144*R^2), c3/(72*R^2),
         c4/(4*R^2)],
    ])


def limit_at_infinity(ff):
    ff = Fn(ff)
    nu = ff.numerator()
    de = ff.denominator()
    dn = nu.degree()
    dd = de.degree()
    if dn < dd:
        return QQ(0)
    if dn == dd:
        return QQ(nu[dn]) / QQ(de[dd])
    raise AssertionError("balanced entry still diverges at infinity: %s" % ff)


u = 2*n + 3
M = authoritative_matrix(u, Fn(R))
D0 = diagonal_matrix(Fn, [1, n, n^2, n^3])
D1 = diagonal_matrix(Fn, [1, n+1, (n+1)^2, (n+1)^3])
B = D0.inverse() * M * D1 / (n+1)^2
S = matrix(QQ, 4, 4, [
    limit_at_infinity(B[i, j]) for i in range(4) for j in range(4)
])

S_expected = matrix(QQ, [
    [64*R-44, 96*R-54, 48*R-17, 8*R],
    [-8, -12, -6, -1],
    [1/R, -4/R, -6/R, -2/R],
    [2/R^2, (17*R-8)/R^2, 4*(5*R-3)/R^2, (6*R-4)/R^2],
])
assert S == S_expected

Rx.<x> = PolynomialRing(QQ)
Q = (
    R^2*x^4
    - (64*R^3 - 56*R^2 - 4)*x^3
    + (48*R^2 - 262*R + 220)*x^2
    - (12*R - 8)*x
    + 1
)
assert Rx(S.charpoly("x")) == Q/R^2
assert Q.is_irreducible()
assert gcd(Q, Q.derivative()) == 1

# On |x|=1 the cubic term strictly dominates all other terms.  Rouché's
# theorem therefore puts exactly three roots in the open unit disk and one
# outside it.
rouche_margin = (
    (64*R^3-56*R^2-4)
    - (R^2 + (48*R^2-262*R+220) + (12*R-8) + 1)
)
assert rouche_margin == 64*R^3-105*R^2+250*R-217
assert rouche_margin > 0

# First row of R^2 adj(xI-S).  At Q(x)=0 it is a left eigenvector.
w = vector(Rx, [
    10 + (44-7*R)*x + (4+12*R^2)*x^2 + R^2*x^3,
    2*((-23+40*R) + (-108+194*R-28*R^2)*x
      + (-27*R^2+48*R^3)*x^2),
    (-32+71*R) + (-68+198*R-8*R^2)*x
      + (-17*R^2+48*R^3)*x^2,
    2*R*(8 + (17+3*R)*x + 4*R^2*x^2),
])
assert w * (x*identity_matrix(Rx, 4) - S.change_ring(Rx)) == vector(Rx, [Q, 0, 0, 0])

# Since Q is irreducible of degree four and every w_j has degree < 4,
# gcd(Q,w_j)=1 proves w_j(rho) != 0 for every root rho of Q.
assert [gcd(Q, z) for z in w] == [Rx(1)]*4

e1 = vector(QQ, [1, 0, 0, 0])
C = matrix(QQ, 4, 4)
for j in range(4):
    C.set_column(j, S^j * e1)
detC_expected = -4*(27*R-11)*(128*R^2-149*R-43)/R^6
assert C.det() == detC_expected
assert C.det() != 0

print("PASS: exact balanced limit and characteristic quartic")
print("PASS: Rouché separation gives three roots inside |x|<1")
print("PASS: explicit left eigenvector has four nonvanishing coordinates")
print("PASS: limiting e1 cyclic frame is invertible")
print("CONCLUSION (using the certified scalar e1 Birkhoff asymptotic):")
print("  lim_N P_(N,j)/Q_(N,j) is independent of j=1,2,3,4")