#!/usr/bin/env python3
"""Account for every family-C survivor. This is the part that matters.

826 four-point operators had vanishing p-curvature at all eight primes I first
tested. None of them is expected to be a counterexample, and the useful output
is a *reason* for each. Two reasons are decidable cheaply and exactly:

  1. More primes. Vanishing at 8 primes is weak evidence; a real candidate must
     vanish at every prime. I re-test each survivor at 24 primes up to 101.
     CPU is free here -- the metered budget on this machine is model tokens, not
     cycles -- so there is no excuse for a weak filter.

  2. Reducibility. If the operator has a rational-function solution it factors,
     and a rank-1 factor with vanishing p-curvature has algebraic solutions
     (rank 1 is the settled case: the residues are rational). So a reducible
     survivor is accounted for and is not a counterexample. I search for a
     polynomial solution exactly over Q -- Fractions, Gaussian elimination, no
     floating point -- because in this family the local exponent at each finite
     singular point is 0, so a rational solution is a polynomial one.

Whatever survives both filters is what I hand over as genuinely unexplained.
"""
import json
from fractions import Fraction as F
from pcurv import vanishes
from sweep import build_C

MORE = [37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97, 101]


def poly_solution(exps, lam, maxdeg=8):
    """Exact search for a polynomial solution of d y'' + Pn y' + (c0 t + c1) y = 0.

    d  = t(t-1)(t-lam),  Pn = e0 (t-1)(t-lam) + e1 t(t-lam) + el t(t-1).
    Returns the coefficient list of a solution, or None.
    """
    e0, e1, el, c0, c1 = [F(*x) for x in exps]
    L = F(*lam)

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

    def add(a, b):
        n = max(len(a), len(b))
        return [(a[i] if i < len(a) else F(0)) + (b[i] if i < len(b) else F(0)) for i in range(n)]

    d = mul(mul([F(0), F(1)], [F(-1), F(1)]), [-L, F(1)])
    Pn = add(add([x * e0 for x in mul([F(-1), F(1)], [-L, F(1)])],
                 [x * e1 for x in mul([F(0), F(1)], [-L, F(1)])]),
             [x * el for x in mul([F(0), F(1)], [F(-1), F(1)])])
    Q = [c1, c0]

    for D in range(0, maxdeg + 1):
        # unknowns y_0..y_D; build the linear system from the coefficient of each
        # power of t in d y'' + Pn y' + Q y.
        rows = D + 4
        M = [[F(0)] * (D + 1) for _ in range(rows)]
        for k in range(D + 1):
            yk = [F(0)] * (D + 1)
            yk[k] = F(1)
            y1 = [F(i) * yk[i] for i in range(1, len(yk))] or [F(0)]
            y2 = [F(i) * y1[i] for i in range(1, len(y1))] or [F(0)]
            col = add(add(mul(d, y2), mul(Pn, y1)), mul(Q, yk))
            for i, v in enumerate(col):
                if i < rows:
                    M[i][k] = v
        # nontrivial kernel?
        ker = kernel(M, D + 1)
        if ker:
            return ker
    return None


def kernel(M, ncols):
    """One nonzero kernel vector of an exact rational matrix, or None."""
    A = [row[:] for row in M]
    nrows = len(A)
    piv, where = 0, [-1] * ncols
    for c in range(ncols):
        r = next((i for i in range(piv, nrows) if A[i][c] != 0), None)
        if r is None:
            continue
        A[piv], A[r] = A[r], A[piv]
        inv = F(1) / A[piv][c]
        A[piv] = [x * inv for x in A[piv]]
        for i in range(nrows):
            if i != piv and A[i][c] != 0:
                f = A[i][c]
                A[i] = [a - f * b for a, b in zip(A[i], A[piv])]
        where[c] = piv
        piv += 1
    free = [c for c in range(ncols) if where[c] == -1]
    if not free:
        return None
    f0 = free[0]
    v = [F(0)] * ncols
    v[f0] = F(1)
    for c in range(ncols):
        if where[c] != -1:
            v[c] = -A[where[c]][f0]
    return v if any(x != 0 for x in v) else None


def main():
    hits = json.load(open('/tmp/hitsC.json'))
    print(f"{len(hits)} survivors from the 8-prime sweep")
    stage2, dead = [], 0
    for (exps, lam, nt) in hits:
        ok = True
        for p in MORE:
            got = build_C((tuple(map(tuple, exps)), tuple(lam)), p)
            if got is None:
                continue
            B, d = got
            if not vanishes(B, d, p):
                ok = False
                break
        if ok:
            stage2.append((exps, lam))
        else:
            dead += 1
    print(f"killed by primes 37..101: {dead}")
    print(f"still vanishing at up to 24 primes: {len(stage2)}")

    unexplained = []
    for (exps, lam) in stage2:
        sol = poly_solution([tuple(x) for x in exps], tuple(lam))
        if sol is None:
            unexplained.append((exps, lam))
    print(f"of those, REDUCIBLE (polynomial solution found): {len(stage2) - len(unexplained)}")
    print(f"of those, no rational solution up to degree 8:   {len(unexplained)}")
    for (exps, lam) in unexplained[:60]:
        e = " ".join(f"{n}/{d}" for (n, d) in exps)
        print(f"   UNEXPLAINED lam={lam[0]}/{lam[1]}  [e0 e1 el c0 c1] = {e}")
    json.dump(unexplained, open('/tmp/unexplainedC.json', 'w'))


if __name__ == "__main__":
    main()