#!/usr/bin/env python3 """The filter I owed the two survivors: symmetric powers of any order, exactly. sym2.py only ruled out the DIHEDRAL case, and I said so in the write-up rather than rounding the result up. An irreducible order-2 operator with finite projective monodromy is dihedral, tetrahedral, octahedral or icosahedral, and those are detected by a rational solution of sym^2, sym^4, sym^6, sym^12 respectively. Both survivors have every exponent a third, which is exactly where the tetrahedral group lives -- so sym^4 is the test that actually decides them, and leaving it unrun was the weakest sentence in the write-up. CONSTRUCTING sym^n L WITHOUT RATIONAL-FUNCTION ARITHMETIC For y'' = -P y' - Q y, put u = y^n for an arbitrary solution y. Every derivative of u lives in the span of the n+1 monomials M_j = y^(n-j) (y')^j, because d/dt M_j = (n-j) M_{j+1} - j P M_j - j Q M_{j-1}. So represent u^(k) as a vector of n+1 coefficients. The coefficients are rational functions, but they are always of the form (polynomial)/d^k with the SAME d, so I never need a general rational-function type: if V/d^m represents a vector, then its derivative is W/d^(m+1) with W_j = (V_j' d - m V_j d') - j Pn V_j + (n-j+1) d V_{j-1} - (j+1) Qn V_{j+1} which is division-free -- the same trick that made the p-curvature sweep cheap. Then u, u', ..., u^(n+1) are n+2 vectors in an (n+1)-dimensional space over Q(t), so a dependence sum_k C_k u^(k) = 0 exists, and its C_k are the coefficients of sym^n L. I find it by scaling everything to the common denominator d^(n+1) and solving one exact linear system for polynomial C_k of bounded degree. WHY THE ORDER MATTERS, AND THE TRAP I AM AVOIDING The test "sym^n L has a rational solution" only means what I want it to mean if the solution space of sym^n L *is* the span of the n-th powers -- i.e. if the order is exactly n+1. If u,...,u^(n) were already dependent, the operator I find is smaller, its solution space is a proper subspace, and a "no rational solution" answer would be about the wrong object. So I check for a dependence among the first n+1 derivatives FIRST and report degeneracy explicitly instead of quietly returning None. This is the same failure I already made once today -- answering a nearby question and having it look like the real one -- and I would rather spend twenty lines than make it twice. """ import json from fractions import Fraction as F from classC import kernel from classC2 import pmul, padd, pscal, pder, parts def binom(n, k): c = 1 for i in range(k): c = c * (n - i) // (i + 1) return c def deriv_vec(V, m, n, d, dd, Pn, Qn): """d/dt of the vector V/d^m in the monomial basis, returned over d^(m+1).""" W = [] for j in range(n + 1): t = padd(pmul(pder(V[j]), d), pscal(pmul(V[j], dd), F(-m)), pscal(pmul(Pn, V[j]), F(-j))) if j >= 1: t = padd(t, pscal(pmul(d, V[j - 1]), F(n - j + 1))) if j + 1 <= n: t = padd(t, pscal(pmul(Qn, V[j + 1]), F(-(j + 1)))) W.append(t) return W def sym_power_operator(exps, lam, n, maxcdeg=40): """Coefficients [C_0..C_{n+1}] of sym^n L, as polynomials over Q. Returns (coeffs, order) or (None, reason). """ d, Pn, Qn, (e0, e1, el), L, (t, t1, tl) = parts(exps, lam) dd = pder(d) N = n + 1 # u = y^n is the vector e_0 over d^0 vecs = [] V = [[F(1)] if j == 0 else [F(0)] for j in range(n + 1)] vecs.append(V) for k in range(N): V = deriv_vec(V, k, n, d, dd, Pn, Qn) vecs.append(V) # scale u^(k) = vecs[k]/d^k up to the common denominator d^N Z = [] for k, V in enumerate(vecs): f = [F(1)] for _ in range(N - k): f = pmul(f, d) Z.append([pmul(c, f) for c in V]) def dependence(cols, maxdeg): """polynomial c_k, not all zero, with sum_k c_k Z[cols[k]] == 0.""" for Dc in range(0, maxdeg + 1): ncols = len(cols) * (Dc + 1) width = max(len(c) for k in cols for c in Z[k]) + Dc + 2 rows = (n + 1) * width M = [[F(0)] * ncols for _ in range(rows)] col = 0 for ci, k in enumerate(cols): for e in range(Dc + 1): for j in range(n + 1): # coefficient polynomial t^e * Z[k][j], placed in block j for i, v in enumerate(Z[k][j]): M[j * width + i + e][col] = v col += 1 v = kernel(M, ncols) # returns ONE vector, not a list of them if v: out = [] for ci in range(len(cols)): out.append([v[ci * (Dc + 1) + e] for e in range(Dc + 1)]) return out return None # trap check: is there already a dependence among u..u^(n)? low = dependence(list(range(n + 1)), min(maxcdeg, 25)) if low: return None, "degenerate: n-th powers span fewer than %d dimensions" % (n + 1) full = dependence(list(range(N + 1)), maxcdeg) if not full: return None, "no polynomial dependence found up to degree %d" % maxcdeg return full, N def rational_solution(coeffs, exps, lam, n, maxdeg=10): """Is there a rational solution of sum_k coeffs[k] u^(k) = 0? A rational solution is prod_s (t-s)^{a_s} * w with w polynomial and a_s a sum of n local exponents at s, i.e. a_s in {k (1 - e_s) : k = 0..n}. """ d, Pn, Qn, (e0, e1, el), L, (t, t1, tl) = parts(exps, lam) dd = pder(d) N = len(coeffs) - 1 opts = [[k * (1 - e) for k in range(n + 1)] for e in (e0, e1, el)] for a0 in opts[0]: for a1 in opts[1]: for al in opts[2]: Rn = padd(pscal(pmul(t1, tl), a0), pscal(pmul(t, tl), a1), pscal(pmul(t, t1), al)) # S_0 = 1, S_{j+1} = S_j' d - j S_j d' + S_j Rn (E^(j)/E = S_j/d^j) S = [[F(1)]] for j in range(N): S.append(padd(pmul(pder(S[j]), d), pscal(pmul(S[j], dd), F(-j)), pmul(S[j], Rn))) dpow = [[F(1)]] for _ in range(N + 1): dpow.append(pmul(dpow[-1], d)) for D in range(0, maxdeg + 1): ncols = D + 1 rows = 2 * (max(len(c) for c in coeffs) + len(dpow[N]) + D + 4) M = [[F(0)] * ncols for _ in range(rows)] for k in range(ncols): wk = [F(0)] * ncols wk[k] = F(1) wd = [wk] for _ in range(N): wd.append(pder(wd[-1])) col = [] for kk in range(N + 1): # U_kk = sum_i C(kk,i) w^(i) S_{kk-i} d^i ; term is C_kk U_kk d^(N-kk) U = [] for i in range(kk + 1): U = padd(U or [F(0)], pscal(pmul(pmul(wd[i], S[kk - i]), dpow[i]), F(binom(kk, i)))) col = padd(col or [F(0)], pmul(pmul(coeffs[kk], U), dpow[N - kk])) for i, v in enumerate(col): if i < rows: M[i][k] = v if kernel(M, ncols): return (a0, a1, al), D return None def check(exps, lam, n, maxdeg=8): coeffs, info = sym_power_operator(exps, lam, n) if coeffs is None: return "sym^%d: could not build the operator (%s)" % (n, info) r = rational_solution(coeffs, exps, lam, n, maxdeg=maxdeg) if r: (a, D) = r return ("sym^%d HAS a rational solution (exponents %s, poly degree <= %d)" " -> finite projective monodromy, solutions algebraic, NOT a counterexample" % (n, tuple(str(x) for x in a), D)) return "sym^%d: no rational solution up to degree %d" % (n, maxdeg) def selftest(): """A control, because a filter that never fires is indistinguishable from a filter that is broken. y'' = ((a+b-1)/t) y' - (ab/t^2) y has solutions t^a and t^b, so y1^n is t^(na) -- rational when na is an integer. In the 4-point family I can force the same thing: exponents that make a power solution rational must be caught by sym^n for the appropriate n.""" print("control: an operator with a rational solution must be caught at sym^1") # e0=e1=el=0, c0=c1=0 -> y'' = 0, solutions 1 and t, both rational. exps = ((0, 1), (0, 1), (0, 1), (0, 1), (0, 1)) for n in (1, 2): print(" n=%d: %s" % (n, check(exps, (2, 1), n, maxdeg=4))) def crosscheck(): """The validation that actually matters, and it needs no outside knowledge. sym2.py builds the symmetric SQUARE from the classical hand-derived formula u''' + 3P u'' + (2P^2 + P' + 4Q) u' + (4PQ + 2Q') u = 0. symn.py builds it from a completely different route: the monomial-basis recursion plus one linear solve, with no formula in it at all. If the two agree as operators -- i.e. their coefficient vectors are proportional -- then the generic machinery is right at n = 2, and n = 4 runs the identical code path with a different n. Without this, a "no rational solution at sym^4" result would be worthless: a broken construction says "no" to everything, and silence is exactly what a negative result looks like. """ from sym2 import ratdiv_free_sym2 # noqa: F401 (kept for provenance) ok = True cases = [(((0, 1), (0, 1), (0, 1), (1, 1), (2, 1)), (2, 1)), (((-2, 3), (-1, 3), (-2, 3), (-1, 1), (1, 1)), (2, 1)), (((1, 2), (2, 3), (-1, 3), (2, 1), (-3, 1)), (3, 1)), (((-1, 3), (-2, 3), (-2, 3), (-1, 1), (0, 1)), (-1, 1))] for exps, lam in cases: d, Pn, Qn, (e0, e1, el), L, _ = parts(exps, lam) dd = pder(d) d2, d3 = pmul(d, d), pmul(pmul(d, d), d) hand = [ 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))), 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))), pscal(pmul(d2, Pn), F(3)), d3, ] gen, info = sym_power_operator(exps, lam, 2) if gen is None: print(" FAIL generic construction failed: %s" % info) ok = False continue # proportional as polynomial vectors: gen[k]*hand[3] == hand[k]*gen[3] good = all(padd(pmul(gen[k], hand[3]), pscal(pmul(hand[k], gen[3]), F(-1))) == [] or all(x == 0 for x in padd(pmul(gen[k], hand[3]), pscal(pmul(hand[k], gen[3]), F(-1)))) for k in range(4)) ok = ok and good e = " ".join(f"{a}/{b}" for (a, b) in exps) print((" pass " if good else " FAIL ") + f"generic sym^2 == hand-derived sym^2 [{e}] lam={lam[0]}/{lam[1]}") print("\npositive control at order 4: e0=3/4, rest 0 -> solutions 1 and t^(1/4),") print("so y^4 = t is rational and sym^4 MUST find a rational solution.") r = check(((3, 4), (0, 1), (0, 1), (0, 1), (0, 1)), (2, 1), 4, maxdeg=4) print(" " + r) ok = ok and "HAS a rational solution" in r print("\nCROSSCHECK " + ("OK" if ok else "FAILED")) return ok if __name__ == "__main__": import sys if "--selftest" in sys.argv: selftest() raise SystemExit(0) if "--crosscheck" in sys.argv: raise SystemExit(0 if crosscheck() else 1) still = json.load(open('/tmp/stillC.json')) orders = [int(x) for x in sys.argv[1:] if x.isdigit()] or [4] print("%d operators survived reducibility, 24 primes and sym^2" % len(still)) for (exps, lam) in still: e = " ".join(f"{a}/{b}" for (a, b) in exps) print("\nlam=%s/%s [%s]" % (lam[0], lam[1], e)) for n in orders: print(" " + check([tuple(x) for x in exps], tuple(lam), n))