#!/usr/bin/env python3
"""Second pass: the reducibility test I ran first was too weak, and here is the fix.

classC.py searched for a *polynomial* solution and called 71 operators
unexplained when it found none. That test is wrong in a way I want on the
record, because it is the same species of error as the wallet-balance one: it
answers a nearby question and it looks like it answered the real one.

An order-2 operator is reducible exactly when the Riccati equation
v' + v^2 + P v + Q = 0 has a RATIONAL solution v = y'/y -- not when y itself is
rational. So an operator with solution (t-1)^(1/2) is reducible, is not a
counterexample to anything, and has no polynomial solution at all. My first
filter would have flagged it as mysterious.

The right finite test, for a Fuchsian operator: a rational v has simple poles
only at the singular points, with residue equal to a local exponent there, plus
negative-integer residues at ordinary points where y has zeros. In this family
the local exponents at each finite singular point s are {0, 1 - e_s} -- Q has
only a simple pole there, so the indicial equation is rho(rho-1) + e_s rho = 0.
So every reducible case has the form

    y = t^m0 (t-1)^m1 (t-lam)^mlam * w(t),   m_s in {0, 1-e_s},   w polynomial

and there are exactly eight choices of exponent vector to try. Substituting and
clearing denominators gives, for R = sum m_s/(t-s) = Rn/d,

    d^2 w'' + d (Pn + 2 Rn) w' + (Qn d + Pn Rn + Rn' d - Rn d' + Rn^2) w = 0

which is a linear system in the coefficients of w, solved exactly over Q.
"""
import json
from fractions import Fraction as F
from classC import kernel


def pmul(a, b):
    out = [F(0)] * (len(a) + len(b) - 1)
    for i, x in enumerate(a):
        if x:
            for j, y in enumerate(b):
                if y:
                    out[i + j] += x * y
    return out


def padd(*ps):
    n = max(len(p) for p in ps)
    return [sum((p[i] if i < len(p) else F(0)) for p in ps) for i in range(n)]


def pscal(a, c):
    return [x * c for x in a]


def pder(a):
    return [F(i) * a[i] for i in range(1, len(a))] or [F(0)]


def parts(exps, lam):
    e0, e1, el, c0, c1 = [F(*x) for x in exps]
    L = F(*lam)
    t, t1, tl = [F(0), F(1)], [F(-1), F(1)], [-L, F(1)]
    d = pmul(pmul(t, t1), tl)
    Pn = padd(pscal(pmul(t1, tl), e0), pscal(pmul(t, tl), e1), pscal(pmul(t, t1), el))
    Qn = [c1, c0]
    return d, Pn, Qn, (e0, e1, el), L, (t, t1, tl)


def reducible(exps, lam, maxdeg=10):
    """(exponent vector, degree of w) for a rational Riccati solution, or None."""
    d, Pn, Qn, (e0, e1, el), L, (t, t1, tl) = parts(exps, lam)
    dd = pder(d)
    for m0 in (F(0), 1 - e0):
        for m1 in (F(0), 1 - e1):
            for ml in (F(0), 1 - el):
                Rn = padd(pscal(pmul(t1, tl), m0), pscal(pmul(t, tl), m1),
                          pscal(pmul(t, t1), ml))
                A2 = pmul(d, d)
                A1 = pmul(d, padd(Pn, pscal(Rn, F(2))))
                A0 = padd(pmul(Qn, d), pmul(Pn, Rn), pmul(pder(Rn), d),
                          pscal(pmul(Rn, dd), F(-1)), pmul(Rn, Rn))
                for D in range(0, maxdeg + 1):
                    rows = D + len(A2) + 2
                    M = [[F(0)] * (D + 1) for _ in range(rows)]
                    for k in range(D + 1):
                        wk = [F(0)] * (D + 1)
                        wk[k] = F(1)
                        w1, w2 = pder(wk), pder(pder(wk))
                        col = padd(pmul(A2, w2), pmul(A1, w1), pmul(A0, wk))
                        for i, v in enumerate(col):
                            if i < rows:
                                M[i][k] = v
                    if kernel(M, D + 1):
                        return ((m0, m1, ml), D)
    return None


def main():
    unexp = json.load(open('/tmp/unexplainedC.json'))
    print(f"{len(unexp)} operators my polynomial-only test could not explain")
    still, explained = [], []
    for (exps, lam) in unexp:
        r = reducible([tuple(x) for x in exps], tuple(lam))
        if r:
            explained.append((exps, lam, r))
        else:
            still.append((exps, lam))
    print(f"reducible after allowing the exponent shift: {len(explained)}")
    for (exps, lam, (m, D)) in explained[:8]:
        e = " ".join(f"{n}/{d}" for (n, d) in exps)
        print(f"   lam={lam[0]}/{lam[1]} [{e}]  y = t^{m[0]} (t-1)^{m[1]} (t-lam)^{m[2]} * poly(deg<={D})")
    print(f"STILL UNEXPLAINED: {len(still)}")
    for (exps, lam) in still:
        e = " ".join(f"{n}/{d}" for (n, d) in exps)
        print(f"   lam={lam[0]}/{lam[1]}  [e0 e1 el c0 c1] = {e}")
    json.dump(still, open('/tmp/stillC.json', 'w'))


if __name__ == "__main__":
    main()