#!/usr/bin/env python3
"""Last filter: the symmetric square, for the two operators that survived everything.

Both survivors have every exponent parameter a third. That is the fingerprint of
finite monodromy, and an irreducible order-2 operator with finite PROJECTIVE
monodromy has a symmetric power with a rational solution -- dihedral shows up
already in the symmetric square. If sym^2 has a rational solution and the
exponents are rational, the solutions are algebraic and the operator is not a
counterexample; it is the four-singular-point analogue of a Schwarz-list case.

For y'' + P y' + Q y = 0 the products y_i y_j satisfy the order-3 equation

    u''' + 3P u'' + (2P^2 + P' + 4Q) u' + (4PQ + 2Q') u = 0

and a rational solution must be prod_s (t-s)^{n_s} times a polynomial, where n_s
is a sum of two local exponents at s, so n_s in {0, 1-e_s, 2(1-e_s)}. Twenty-seven
exponent vectors, one exact linear solve each.
"""
import json
from fractions import Fraction as F
from classC import kernel
from classC2 import pmul, padd, pscal, pder, parts


def ratdiv_free_sym2(exps, lam, maxdeg=12):
    d, Pn, Qn, (e0, e1, el), L, (t, t1, tl) = parts(exps, lam)
    # P = Pn/d, Q = Qn/d. Clear denominators: multiply the order-3 equation by d^3.
    dd = pder(d)
    d2, d3 = pmul(d, d), pmul(pmul(d, d), d)
    # coefficients of u''', u'', u', u after multiplying by d^3
    C3 = d3
    C2 = pscal(pmul(d2, Pn), F(3))
    # (2P^2 + P' + 4Q) d^3 = 2 Pn^2 d + (Pn' d - Pn d') d + 4 Qn d^2
    C1 = padd(pscal(pmul(pmul(Pn, Pn), d), F(2)),
              pmul(padd(pmul(pder(Pn), d), pscal(pmul(Pn, dd), F(-1))), d),
              pscal(pmul(Qn, d2), F(4)))
    # (4PQ + 2Q') d^3 = 4 Pn Qn d + 2 (Qn' d - Qn d') d
    C0 = padd(pscal(pmul(pmul(Pn, Qn), d), F(4)),
              pscal(pmul(padd(pmul(pder(Qn), d), pscal(pmul(Qn, dd), F(-1))), d), F(2)))

    for n0 in (F(0), 1 - e0, 2 * (1 - e0)):
        for n1 in (F(0), 1 - e1, 2 * (1 - e1)):
            for nl in (F(0), 1 - el, 2 * (1 - el)):
                # u = t^n0 (t-1)^n1 (t-lam)^nl * w  ->  shift by R = sum n_s/(t-s)
                Rn = padd(pscal(pmul(t1, tl), n0), pscal(pmul(t, tl), n1),
                          pscal(pmul(t, t1), nl))
                # Substituting u = E w with E'/E = Rn/d into the order-3 operator
                # and clearing d^3 gives coefficients built from Rn. Rather than
                # expand by hand, build the operator's action on w numerically in
                # the polynomial ring: (E w)^(k) / E expressed via Rn/d.
                # s1 = Rn/d, s2 = (Rn/d)' + (Rn/d)^2, s3 = ...
                # Work with numerators over d, d^2, d^3 respectively.
                a1 = Rn                                              # /d
                a2 = padd(pmul(pder(Rn), d), pscal(pmul(Rn, dd), F(-1)), pmul(Rn, Rn))  # /d^2
                # a3 = a2' d - 2 a2 d' + a1 a2   (numerator over d^3)
                a3 = padd(pmul(pder(a2), d), pscal(pmul(a2, dd), F(-2)), pmul(a1, a2))
                for D in range(0, maxdeg + 1):
                    ncols = D + 1
                    rows = D + len(C3) + 4
                    M = [[F(0)] * ncols for _ in range(rows)]
                    for k in range(ncols):
                        wk = [F(0)] * ncols
                        wk[k] = F(1)
                        w1, w2, w3 = pder(wk), pder(pder(wk)), pder(pder(pder(wk)))
                        # (E w)''' / E * d^3 = w''' d^3 + 3 w'' a1 d^2 + 3 w' a2 d + w a3
                        u3 = padd(pmul(w3, d3), pscal(pmul(pmul(w2, a1), d2), F(3)),
                                  pscal(pmul(pmul(w1, a2), d), F(3)), pmul(wk, a3))
                        u2 = padd(pmul(w2, d2), pscal(pmul(pmul(w1, a1), d), F(2)), pmul(wk, a2))
                        u1 = padd(pmul(w1, d), pmul(wk, a1))
                        # C3 u3/d^3 + C2 u2/d^2 + C1 u1/d + C0 u = 0, times d^3/d^3...
                        col = padd(pmul(C3, u3), pmul(pmul(C2, u2), d),
                                   pmul(pmul(C1, u1), d2), pmul(pmul(C0, wk), d3))
                        for i, v in enumerate(col):
                            if i < rows:
                                M[i][k] = v
                    if kernel(M, ncols):
                        return (n0, n1, nl), D
    return None


if __name__ == "__main__":
    still = json.load(open('/tmp/stillC.json'))
    print(f"{len(still)} operators survived reducibility and 24 primes")
    for (exps, lam) in still:
        e = " ".join(f"{n}/{d}" for (n, d) in exps)
        r = ratdiv_free_sym2([tuple(x) for x in exps], tuple(lam))
        print(f"lam={lam[0]}/{lam[1]}  [{e}]")
        if r:
            (n, D) = r
            print(f"   sym^2 HAS a rational solution: exponents {n}, poly degree <= {D}")
            print("   -> projective monodromy dihedral/finite, solutions algebraic, NOT a counterexample")
        else:
            print("   sym^2 has NO rational solution up to degree 12")
            print("   -> not explained by any filter I have. Handing this one over.")