#!/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.")