#!/usr/bin/env python3 """ certificate_irreducibility.py -- Lemma "Irreducibility over the base field" (the first of the four lemmas proving the no-go proposition): L5 has no right factor of order 1, 2, 3, or 4 over Q(xi). L5 is the order-five Picard-Fuchs operator of the orthorhombic Watson integral at the physical rates, (A,B,C) = (alpha^2,beta^2,gamma^2) = (2,3,5), acting on w W_S as a function of xi = 1/w^2: the operator listed in the appendix, stored in theta = xi d/dxi form as data/L5_theta.txt. Every identity below is exact integer or rational arithmetic; the one modular step is a rank bound that is rigorous in the direction used (rank mod p <= rank over Q), and it is used only to certify that a linear system has NO solution. No floating-point number enters anywhere. STAGES (each prints PASS/FAIL lines; the certificate is fail-closed: any FAIL, any completeness flag, or any exception gives "CERTIFICATE FAIL" and a nonzero exit status). Stage 0 The operator. data/L5_theta.txt is specialized to (A,B,C)=(2,3,5) and converted from theta form to the D = d/dxi form over Z[xi] (theta^n = sum_k S(n,k) xi^k D^k with Stirling numbers of the second kind; denominators and content cleared), and the result must equal the stored coefficient lists [c_0,...,c_5] of data/L5_dform.json that Stages 1-3 consume. Its leading coefficient is c_5 = xi^5 (216 xi^2 + 20 xi + 3)(576 xi^4 - 960 xi^3 + 352 xi^2 - 40 xi + 1), the singular locus of the paper (xi = 0, the apparent quadratic, the threshold quartic, infinity). Stage 1 Order one: complete hyperexponential enumeration for L5 itself. A right factor of order one over Q(xi) is D - u'/u with u a hyperexponential solution of L5, u = prod_F F^{e_F} P(xi) with F running over the irreducible factors of c_5 over Q, e_F the exponent at the roots of F (a root of the indicial polynomial there, the same at every root of F: the exponents are the residues of u'/u in Q(xi), Galois conjugate over the roots of F, and rational at a rational point; a rational exponent is enumerated, an irrational one at a factor of degree > 1 raises a completeness flag), nonnegative-integer exponents absorbed into the polynomial P, and deg P = -lam - sum_F e_F deg F pinned by a root lam of the indicial polynomial at infinity (u ~ xi^{-lam}). The certificate lists every exponent option at every factor and every lam, and checks (i) that the enumeration engine finds no combination with a nonnegative integer degree, (ii) directly that every finite exponent option is >= 0 and every lam is > 0, so that every combination forces deg P < 0 (the sentence in the paper), and (iii) that the roots at infinity are those of the printed P_6, {1/2, 1, 3/2, 5/2, 7/2}. Completeness conditions are FAIL conditions: at every factor of c_5 the indicial polynomial must have constant coefficients and all of its roots found by exact factorization over Q, and at infinity every root must be rational and found. PASS = no hyperexponential solution, hence no order-one right factor. Stage 2 The exterior-power annihilators are exact. data/exterior_powers.json stores operators M_2, M_3, M_4 over Z[xi] (orders 10, 10, 5; coefficient degrees 76, 90, 29; integer coefficients of at most 233, 277, 90 bits). For each r the certificate proves, in exact integer-polynomial arithmetic with no series truncation, that M_r annihilates the r x r Wronskian of EVERY r-tuple of solutions of L5. With m the vector of r-minors of the derivative matrix (index sets I in {0,...,4}, |I| = r), m' = B m is the connection induced by the companion matrix of L5, B = Bnum/c_5; the derivative tower (d/dxi)^k m = B_k m, B_k = N_k / c_5^k, obeys N_{k+1} = N_k' c_5 - k N_k c_5' + N_k Bnum in integer-polynomial matrices, and M_r = sum_k ct_k D^k annihilates the minor m_{(0,...,r-1)} of every solution tuple iff row (0,...,r-1) of sum_k ct_k N_k c_5^{n-k} is identically zero in Z[xi]. PASS = that row is the zero vector (residual degree -1) for r = 2, 3, 4. How the M_r were found (series fits modulo primes, then rational reconstruction) plays no role in the proof: only the exact identity certified here is used, and minimality of M_r is not needed for Stage 3. Stage 3 Orders two, three, four: Beke's criterion on M_r. If L5 = Q R with R of order r over Q(xi), the Wronskian of the solution space of R is hyperexponential (its logarithmic derivative is minus the ratio of the two top coefficients of R, in Q(xi)) and it is an r-wedge of solutions of L5, hence a hyperexponential solution of M_r. The Stage 1 enumeration is therefore run on M_2, M_3, M_4: factor the leading coefficient over Q, take the exponent options from the indicial roots at every factor, pin deg P by the indicial roots at infinity, and for each admissible combination set up the finite linear system over Q for the coefficients of P. NO solution exists iff that matrix has full column rank; full column rank is certified modulo a 61-bit prime (a second prime is consulted only if the first is rank-deficient; rank mod p <= rank over Q, so full rank mod p proves full rank over Q). PASS = every combination is NOSOL, no candidate, no completeness flag, and the same root census as in Stage 1 at every factor. For r = 2 and 3 the annihilators have order ten; for r = 4 the exact order-five integer-coefficient annihilator has two admissible combinations, both excluded by full column rank. WHAT PASS ESTABLISHES. L5 has no right factor of order 1, 2, 3 or 4 in Q(xi): the lemma as printed. Since any proper factorization of an order-five operator has a right factor of order 1 to 4, L5 is irreducible over Q(xi). Descent to finite extensions of Q(xi), the local structure at xi = 0, and the exclusion of symmetric powers are the other three lemmas of the proposition and are certified elsewhere (frobenius_L5.py for the local structure); they are not part of this file. Inputs: data/L5_theta.txt, data/L5_dform.json, data/exterior_powers.json (SHA-256 of each is recorded in the JSON verdict). Outputs: certificate_irreducibility.log (this report) and certificate_irreducibility_result.json (verdict, per-stage data). Expected: "SUMMARY: 24 PASS, 0 FAIL" then "CERTIFICATE PASS (fail-closed)"; about 15 s (the M_3 sweep dominates). Requires: Python >= 3.10, sympy, python-flint (fmpz_poly, fmpq_poly, nmod_mat). """ import sys, json, time, hashlib from fractions import Fraction as Fr from itertools import combinations, product as iproduct import sympy as sp import flint # ========================================================================= # Stage 0 support: the operator, theta form -> D form over Z[xi]. # ========================================================================= # --- theta_symbols: begin --- x, th = sp.symbols('x theta') A, B, C = sp.symbols('A B C') # --- theta_symbols: end --- def dform_from_theta(): """Return [c_0,...,c_5] as integer coefficient lists (low -> high degree) of L5 in D-form over Z[xi] at (A,B,C) = (2,3,5), built from data/L5_theta.txt.""" # --- theta_to_D: begin --- # parse L5_theta.txt ops = {} for line in open('data/L5_theta.txt'): line = line.strip() if not line or line.startswith('#'): continue lhs, rhs = line.split(':', 1) p = int(lhs.replace('x^', '')) ops[p] = sp.sympify(rhs, locals={'theta': th, 'A': A, 'B': B, 'C': C}) subs = {A: 2, B: 3, C: 5} # theta^k -> sum_j S2(k,j) x^j D^j D = sp.symbols('D') # placeholder, track as list maxord = 5 # build c_i(x): operator = sum_p x^p P_p(theta); theta^n = sum_k S(n,k) x^k D^k from sympy.functions.combinatorial.numbers import stirling coef = [sp.Integer(0)] * (maxord + 1) for p, P in ops.items(): Pn = sp.Poly(sp.expand(P.subs(subs)), th) for (n,), a in Pn.terms(): for k in range(n + 1): coef[k] += a * stirling(n, k, kind=2) * x**(p + k) # n=0 handled: stirling(0,0)=1 coef = [sp.expand(c) for c in coef] # clear denominators den = sp.lcm([sp.nsimplify(sp.denom(sp.together(t))) for c in coef for t in sp.Add.make_args(c)] + [1]) coef = [sp.expand(c * den) for c in coef] g = sp.gcd(list(sp.Poly(sum(coef[i]*sp.Symbol(f'D{i}') for i in range(6)), x, *[sp.Symbol(f'D{i}') for i in range(6)]).coeffs())) coef = [sp.expand(c / g) for c in coef] cpolys = [sp.Poly(c, x) for c in coef] print("degrees:", [p.degree() if not p.is_zero else -1 for p in cpolys]) print("lc(L5) =", sp.factor(coef[5])) # integer lists low->high il = [] for p in cpolys: d = p.degree() if not p.is_zero else 0 il.append([int(p.coeff_monomial(x**j)) for j in range(d + 1)]) # --- theta_to_D: end --- return il # ========================================================================= # Stage 2 engine: the exact wedge-annihilation identity (companion matrix of # L5, induced connection on r-minors, integer-polynomial derivative tower). # ========================================================================= L5 = json.load(open('data/L5_dform.json')) EXT = json.load(open('data/exterior_powers.json')) # --- companion_setup: begin --- Z = flint.fmpz_poly lc = Z(L5[5]) lcd = lc.derivative() # companion numerator: (y^(i))' = sum_j Cn[i][j]/lc * y^(j) Cn = [[Z([0])] * 5 for _ in range(5)] for i in range(4): Cn[i][i + 1] = lc for j in range(5): Cn[4][j] = Z([-c for c in L5[j]]) # --- companion_setup: end --- # --- wedge_identity: begin --- def wedge_Bnum(r): basis = list(combinations(range(5), r)) idx = {I: i for i, I in enumerate(basis)} nb = len(basis) Bn = [[Z([0])] * nb for _ in range(nb)] for I in basis: for kpos, ik in enumerate(I): for j in range(5): if Cn[ik][j].degree() < 0 or Cn[ik][j] == Z([0]): continue newI = list(I) newI[kpos] = j if len(set(newI)) < r: continue srt = tuple(sorted(newI)) # permutation sign to sort newI sign = 1 arr = list(newI) for a in range(len(arr)): for b in range(a + 1, len(arr)): if arr[a] > arr[b]: arr[a], arr[b] = arr[b], arr[a] sign = -sign Bn[idx[I]][idx[srt]] += Cn[ik][j] if sign == 1 else -Cn[ik][j] return basis, Bn def prove(r): op = EXT[str(r)] n = len(op) - 1 ct = [Z(ci) for ci in op] basis, Bn = wedge_Bnum(r) nb = len(basis) I0 = tuple(range(r)) i0 = basis.index(I0) # row-vector tower: v_k = row i0 of N_k (N_0 = Id) v = [Z([0])] * nb v[i0] = Z([1]) # accumulate row of sum_k ct_k N_k lc^{n-k} lcpow = [Z([1])] for _ in range(n): lcpow.append(lcpow[-1] * lc) acc = [ct[0] * lcpow[n] * vv for vv in v] for k in range(1, n + 1): # v <- v' * lc - (k-1) * v * lc' + v * Bnum (row update: (N v as row) row_{k} = row_{k-1}'*lc - (k-1)row_{k-1}lc' + row_{k-1}*Bnum) nv = [Z([0])] * nb for j in range(nb): nv[j] = v[j].derivative() * lc - v[j] * lcd * (k - 1) for j in range(nb): if v[j] == Z([0]): continue for l in range(nb): if Bn[j][l] != Z([0]): nv[l] += v[j] * Bn[j][l] v = nv for j in range(nb): acc[j] += ct[k] * lcpow[n - k] * v[j] zero = all(a == Z([0]) for a in acc) mx = max((a.degree() for a in acc), default=-1) print(f"r={r}: order {n}, wedge dim {nb}, row-identity zero: {zero} (max residual deg {mx})") return zero # --- wedge_identity: end --- # ========================================================================= # Stages 1 and 3 engine: complete hyperexponential enumeration (indicial # polynomials at every leading-coefficient factor and at infinity, exponent # options, degree pinning, and the full-column-rank exclusion per combination). # ========================================================================= # --- sweep_constants: begin --- P1 = (1 << 61) - 1 P2 = 2305843009213693967 xs = sp.symbols('x') # --- sweep_constants: end --- # --- sweep_engine: begin --- def il_to_sympy(ci): return sp.Poly(list(reversed(ci)), xs) def sympy_to_fmpq(poly): d = poly.degree() if not poly.is_zero else 0 return flint.fmpq_poly([sp.Rational(poly.coeff_monomial(xs**j)) and Fr(int(sp.nsimplify(poly.coeff_monomial(xs**j)).p), int(sp.nsimplify(poly.coeff_monomial(xs**j)).q)) or Fr(0) for j in range(d + 1)]) def fq(fr_list): return flint.fmpq_poly([flint.fmpq(f.numerator, f.denominator) for f in fr_list]) def fq_from_ints(il): return flint.fmpq_poly(il) def indicial_at_factor(cs_sym, Fpoly): """Indicial polynomial of M at roots s of irreducible Fpoly. Returns sympy Poly in theta with coefficients in Q[s]/F (as sympy exprs in s), plus the multiplicity vector v_i.""" s = sp.symbols('s') n = len(cs_sym) - 1 v = [] g_mod = [] for ci in cs_sym: vi = 0 rem = ci while True: q, r = sp.div(rem, Fpoly, xs) if r.is_zero and not rem.is_zero: vi += 1 rem = q else: break v.append(vi) g_mod.append(sp.rem(rem, Fpoly, xs)) # g_i mod F k = max(i - v[i] for i in range(n + 1) if not cs_sym[i].is_zero) th = sp.symbols('theta') Fp_der = sp.diff(Fpoly.as_expr(), xs) ind = sp.Integer(0) for i in range(n + 1): if cs_sym[i].is_zero or i - v[i] != k: continue fall = sp.prod([th - j for j in range(i)]) term = g_mod[i].as_expr() * sp.rem(sp.Poly(Fp_der**v[i], xs), Fpoly, xs).as_expr() * fall ind += term # reduce coefficients mod F indp = sp.Poly(sp.expand(ind), th) red = [] for c in indp.all_coeffs(): red.append(sp.rem(sp.Poly(c, xs), Fpoly, xs).as_expr()) return sp.Poly(red, th), v, k def exponent_options(indpoly_theta, Fpoly): """Roots of indicial poly that are rational constants; classify. Returns (options, complete_flag). options = list of Fraction exponents to try (0 always included; nonneg integers dropped as absorbable).""" th = indpoly_theta.gens[0] # coefficients must be constants (degree 0 in x) after normalization coeffs = indpoly_theta.all_coeffs() frees = set() for c in coeffs: frees |= sp.sympify(c).free_symbols complete = True opts = {Fr(0)} if frees - {th}: # non-constant coefficients in Q[s]/F: factor over extension needed # try: normalize monic and re-check lead = coeffs[0] norm = [sp.simplify(sp.cancel(c / lead)) for c in coeffs] frees2 = set() for c in norm: frees2 |= sp.sympify(c).free_symbols if frees2 - {th}: return None, False # caller must handle via extension factoring coeffs = norm ip = sp.Poly(coeffs, th) rts = sp.roots(ip) unaccounted = 0 for rt, mult in rts.items(): if rt.is_rational: f = Fr(int(rt.p), int(rt.q)) if f.denominator == 1 and f >= 0: continue # absorbed into P opts.add(f) else: unaccounted += mult if unaccounted and sum(rts.values()) < ip.degree(): unaccounted += ip.degree() - sum(rts.values()) # irrational roots: can they lie in Q(s)? only if deg F > 1; flag incomplete if unaccounted and Fpoly.degree() > 1: complete = 'IRRATIONAL_ROOTS_NEED_EXTENSION_CHECK' return sorted(opts), complete def indicial_at_infinity(cs_sym): """ind_inf(lam) = sum_{i: deg c_i - i = m*} lc(c_i) * (-1)^i lam(lam+1)...(lam+i-1). Exponent lam means w ~ x^{-lam} at infinity.""" n = len(cs_sym) - 1 degs = [(-10**9 if c.is_zero else c.degree()) for c in cs_sym] mstar = max(degs[i] - i for i in range(n + 1)) lam = sp.symbols('lam') ind = sp.Integer(0) for i in range(n + 1): if degs[i] - i == mstar: rise = sp.prod([lam + j for j in range(i)]) ind += cs_sym[i].LC() * (-1)**i * rise return sp.Poly(sp.expand(ind), lam) def lam_options(indinf): lamr = sp.roots(indinf) opts = [] complete = True tot = 0 for rt, mult in lamr.items(): tot += mult if rt.is_rational: opts.append(Fr(int(rt.p), int(rt.q))) else: complete = 'IRRATIONAL_LAMBDA' if tot < indinf.degree(): complete = 'MISSING_LAMBDA_ROOTS' return opts, complete def build_Gj(cs_fq, Ud, Un, n): """G_j polys (fmpq_poly) with sum_j G_j P^(j) = 0 <=> M(prod F^e * P) = 0. T_{0,0}=1; T_{k+1,j} = T'_{k,j} Ud + T_{k,j-1} Ud - k T_{k,j} Ud' + Un T_{k,j}. G_j = sum_i c_i Ud^{n-i} T_{i,j}.""" Udp = Ud.derivative() T = {(0, 0): flint.fmpq_poly([1])} for k in range(n): for j in range(k + 2): t = flint.fmpq_poly([0]) if (k, j) in T: t += T[(k, j)].derivative() * Ud - T[(k, j)] * Udp * k + Un * T[(k, j)] if (k, j - 1) in T: t += T[(k, j - 1)] * Ud if t.length() or True: T[(k + 1, j)] = t # powers of Ud Upow = [flint.fmpq_poly([1])] for i in range(n): Upow.append(Upow[-1] * Ud) G = [] for j in range(n + 1): g = flint.fmpq_poly([0]) for i in range(j, n + 1): if (i, j) in T: g += cs_fq[i] * Upow[n - i] * T[(i, j)] G.append(g) return G def combo_matrix(G, D): """Rows: coefficients of x^m of sum_j G_j d^j/dx^j x^k, k=0..D. Fraction matrix.""" n = len(G) - 1 maxlen = max(g.length() for g in G) + D + 2 cols = [] for k in range(D + 1): col = [Fr(0)] * maxlen for j in range(min(k, n) + 1): fall = 1 for t in range(j): fall *= (k - t) if fall == 0: continue gj = G[j] for m in range(gj.length()): c = gj[m] idx = m + k - j if idx < maxlen: col[idx] += Fr(int(c.p), int(c.q)) * fall # flint fmpz -> int (Fraction rejects fmpz) cols.append(col) return cols, maxlen def rank_modp(cols, nrows, p): ncols = len(cols) flat = [] # row-major den_scale = [] ints = [] for col in cols: den = 1 for c in col: den = den * c.denominator // sp.igcd(den, c.denominator) ints.append([int(c * den) for c in col]) for m in range(nrows): for k in range(ncols): flat.append(ints[k][m] % p) M = flint.nmod_mat(nrows, ncols, flat, p) return M.rank(), ints def sweep_operator(name, op_il, results): t0 = time.time() n = len(op_il) - 1 cs_sym = [il_to_sympy(c) for c in op_il] lc = cs_sym[n] fl = sp.factor_list(lc.as_expr()) factors = [(sp.Poly(f, xs), m) for f, m in fl[1] if sp.Poly(f, xs).degree() > 0] print(f"[{name}] order {n}, lc factors: {[ (str(f.as_expr())[:40], m) for f,m in factors ]}") # Fuchs check (info): v_i multiplicities opt_sets = [] flags = [] for F, mult in factors: ind, v, k = indicial_at_factor(cs_sym, F) opts, complete = exponent_options(ind, F) if opts is None: print(f" factor {F.as_expr()}: NON-CONSTANT indicial coeffs — needs extension handling") flags.append((str(F.as_expr()), 'NONCONST')) opts, complete = [Fr(0)], 'NONCONST_TRUNCATED' if complete is not True: flags.append((str(F.as_expr()), str(complete))) print(f" factor deg{F.degree()}: v={v} k={k} exponent options {opts} complete={complete}") opt_sets.append((F, opts)) indinf = indicial_at_infinity(cs_sym) lams, lcomplete = lam_options(indinf) if lcomplete is not True: flags.append(('infinity', str(lcomplete))) print(f" infinity: lam options {sorted(lams)} complete={lcomplete}") cs_fq = [fq_from_ints(c) for c in op_il] Ud = flint.fmpq_poly([1]) for F, _ in factors: d = F.degree() Ud = Ud * flint.fmpq_poly([flint.fmpq(sp.Rational(F.coeff_monomial(xs**j)).p, sp.Rational(F.coeff_monomial(xs**j)).q) for j in range(d + 1)]) Fpolys = [flint.fmpq_poly([flint.fmpq(sp.Rational(F.coeff_monomial(xs**j)).p, sp.Rational(F.coeff_monomial(xs**j)).q) for j in range(F.degree() + 1)]) for F, _ in factors] ncombo = 0 candidates = [] for choice in iproduct(*[opts for _, opts in opt_sets]): for lam in lams: D = -lam - sum(e * F.degree() for (F, _), e in zip(opt_sets, choice)) if D.denominator != 1 or D < 0: continue D = int(D) if D > 400: flags.append((str(choice), f'D={D}>400 SKIPPED')) continue ncombo += 1 # Un = sum_F e_F F' * (Ud/F) Un = flint.fmpq_poly([0]) for i, ((F, _), e) in enumerate(zip(opt_sets, choice)): if e == 0: continue cof = Ud // Fpolys[i] Un += Fpolys[i].derivative() * cof * flint.fmpq_poly([flint.fmpq(e.numerator, e.denominator)]) G = build_Gj(cs_fq, Ud, Un, n) cols, nrows = combo_matrix(G, D) r1, ints = rank_modp(cols, nrows, P1) status = 'NOSOL(rank p1)' if r1 < D + 1: r2, _ = rank_modp(cols, nrows, P2) if r2 < D + 1: status = 'CANDIDATE' candidates.append((choice, lam, D)) else: status = 'NOSOL(rank p2)' print(f" combo e={[str(e) for e in choice]} lam={lam} D={D}: {status}") out = {'operator': name, 'order': n, 'ncombos_solved': ncombo, 'candidates': [( [str(e) for e in c], str(l), D) for c, l, D in candidates], 'flags': flags, 'seconds': round(time.time() - t0, 1)} results[name] = out print(f"[{name}] done: {ncombo} combos, {len(candidates)} candidates, flags={flags} ({out['seconds']}s)") return out # --- sweep_engine: end --- # ========================================================================= # Driver: PASS/FAIL bookkeeping, log tee, and the three stages. # ========================================================================= LOGNAME = 'certificate_irreducibility.log' RESNAME = 'certificate_irreducibility_result.json' class _Tee: def __init__(self, *streams): self.streams = streams def write(self, s): for st in self.streams: st.write(s) def flush(self): for st in self.streams: st.flush() npass = nfail = 0 checks = [] def gate(name, cond, detail=''): global npass, nfail ok = bool(cond) if ok: npass += 1 else: nfail += 1 checks.append({'check': name, 'pass': ok, 'detail': detail}) print(('PASS ' if ok else 'FAIL ') + name + ((' -- ' + detail) if detail else '')) return ok def sha256(path): return hashlib.sha256(open(path, 'rb').read()).hexdigest() def indicial_tables(op_il): """Exponent options at every leading-coefficient factor and at infinity, computed with the sweep's own functions (same calls as sweep_operator).""" n = len(op_il) - 1 cs_sym = [il_to_sympy(c) for c in op_il] lc = cs_sym[n] fl = sp.factor_list(lc.as_expr()) factors = [(sp.Poly(f, xs), m) for f, m in fl[1] if sp.Poly(f, xs).degree() > 0] rows = [] for F, mult in factors: ind, v, k = indicial_at_factor(cs_sym, F) opts, complete = exponent_options(ind, F) rows.append((F, mult, v, k, opts, complete, ind)) indinf = indicial_at_infinity(cs_sym) lams, lcomplete = lam_options(indinf) return cs_sym, factors, rows, sorted(lams), lcomplete def root_census(rows): """Completeness guard over and above the sweep's flags: at every irreducible factor F of the leading coefficient the indicial polynomial, normalized to be monic, must have constant coefficients and ALL of its roots found by exact factorization over Q (sum of multiplicities equal to its degree), so that no rational exponent can have been missed and no irrational exponent at a factor of degree > 1 can have gone unflagged.""" ok = True table = [] for F, mult, v, k, opts, complete, ind in rows: th = ind.gens[0] coeffs = ind.all_coeffs() lead = coeffs[0] norm = [sp.simplify(sp.cancel(c / lead)) for c in coeffs] frees = set() for c in norm: frees |= sp.sympify(c).free_symbols if frees - {th}: ok = False table.append({'factor': str(F.as_expr())[:40], 'deg': F.degree(), 'status': 'NONCONSTANT indicial coefficients'}) continue ip = sp.Poly(norm, th) rts = sp.roots(ip) found = sum(rts.values()) allfound = (found == ip.degree()) ok = ok and allfound table.append({'factor': str(F.as_expr())[:40], 'deg': F.degree(), 'indicial_degree': ip.degree(), 'roots_found': found, 'roots': sorted(str(r) for r in rts), 'status': 'all roots found' if allfound else 'MISSING ROOTS'}) return ok, table def main(): t0 = time.time() LOG = open(LOGNAME, 'w') sys.stdout = _Tee(sys.__stdout__, LOG) result = {'lemma': 'L5 has no right factor of order 1, 2, 3, or 4 over Q(xi)', 'rates': '(A,B,C) = (2,3,5)', 'inputs': {}, 'software': {}} for f in ('data/L5_theta.txt', 'data/L5_dform.json', 'data/exterior_powers.json'): result['inputs'][f] = sha256(f) result['software'] = {'python': sys.version.split()[0], 'sympy': sp.__version__, 'python-flint': flint.__version__} print('certificate_irreducibility.py -- L5 at (A,B,C) = (2,3,5): no right factor of order 1-4 over Q(xi)') for f, h in result['inputs'].items(): print(f'input {f}: sha256 {h}') print('software:', result['software']) # ---------------- Stage 0: the operator ---------------- print('\n== Stage 0: operator (theta form of data/L5_theta.txt -> D form over Z[xi]) ==') il = dform_from_theta() L5_stored = json.load(open('data/L5_dform.json')) gate('Stage 0: D-form of data/L5_theta.txt at (2,3,5) == data/L5_dform.json', il == L5_stored) lcf = sp.factor(il_to_sympy(L5_stored[5]).as_expr()) lc_expected = sp.factor(xs**5 * (216 * xs**2 + 20 * xs + 3) * (576 * xs**4 - 960 * xs**3 + 352 * xs**2 - 40 * xs + 1)) gate('Stage 0: leading coefficient c_5 = xi^5 (216 xi^2+20 xi+3)(576 xi^4-960 xi^3+352 xi^2-40 xi+1)', sp.expand(lcf - lc_expected) == 0) gate('Stage 0: L5_dform.json has order 5 and integer coefficient lists', len(L5_stored) == 6 and all(isinstance(c, int) for ci in L5_stored for c in ci)) gate('Stage 0: certificate operators agree (L5 of prove_wedge == L5_dform.json)', L5 == L5_stored) result['stage0'] = {'degrees': [len(c) - 1 for c in L5_stored], 'lc_factored': str(lcf)} # ---------------- Stage 1: order one ---------------- print('\n== Stage 1: order one, complete hyperexponential enumeration for L5 ==') results = {} r1 = sweep_operator('L5_r1', L5_stored, results) gate('Stage 1: engine finds no exponent combination with a nonnegative integer degree (0 combinations, 0 candidates)', r1['ncombos_solved'] == 0 and r1['candidates'] == []) gate('Stage 1: no completeness flag raised by the enumeration', r1['flags'] == []) cs_sym, factors, rows, lams, lcomplete = indicial_tables(L5_stored) print(' exponent options per factor of c_5:') for F, mult, v, k, opts, complete, ind in rows: print(f' factor deg {F.degree()} ({str(F.as_expr())[:40]}): options {[str(e) for e in opts]} complete={complete}') print(f' roots at infinity (u ~ xi^-lam): {[str(l) for l in lams]} complete={lcomplete}') all_e_nonneg = all(e >= 0 for _, _, _, _, opts, _, _ in rows for e in opts) all_lam_pos = all(l > 0 for l in lams) gate('Stage 1: every finite exponent option >= 0 and every lam > 0, so deg P = -lam - sum e_F deg F < 0 for every combination', all_e_nonneg and all_lam_pos and all(c is True for _, _, _, _, _, c, _ in rows) and lcomplete is True) gate('Stage 1: roots at infinity are those of the printed P_6: {1/2, 1, 3/2, 5/2, 7/2}', lams == [Fr(1, 2), Fr(1), Fr(3, 2), Fr(5, 2), Fr(7, 2)]) ok_c, census = root_census(rows) gate('Stage 1: root census -- at every factor of c_5 the indicial polynomial has constant coefficients and all roots found over Q', ok_c) for row in census: print(' ', row) result['stage1'] = {'sweep': r1, 'exponent_options': [{'factor': str(F.as_expr())[:40], 'deg': F.degree(), 'options': [str(e) for e in opts], 'complete': str(c)} for F, mult, v, k, opts, c, ind in rows], 'lam_options': [str(l) for l in lams], 'root_census': census} # ---------------- Stage 2: exact exterior-power annihilators ---------------- print('\n== Stage 2: M_r annihilates every r-wedge of solutions of L5, exactly over Z[xi] ==') expected_order = {'2': 10, '3': 10, '4': 5} result['stage2'] = {} for r in (2, 3, 4): op = EXT[str(r)] n = len(op) - 1 maxbits = max(abs(int(c)).bit_length() for ci in op for c in ci) gate(f'Stage 2: M_{r} has order {expected_order[str(r)]} with integer coefficients', n == expected_order[str(r)] and all(isinstance(c, int) for ci in op for c in ci), f'order {n}, max coefficient degree {max(len(c) - 1 for c in op)}, max {maxbits} bits') ok = prove(r) gate(f'Stage 2: exact wedge annihilation identity for r = {r} (row identity zero in Z[xi])', ok) result['stage2'][f'M_{r}'] = {'order': n, 'max_degree': max(len(c) - 1 for c in op), 'max_bits': maxbits, 'annihilates_all_wedges': bool(ok)} # ---------------- Stage 3: Beke sweeps on M_2, M_3, M_4 ---------------- print('\n== Stage 3: complete hyperexponential enumeration for M_2, M_3, M_4 ==') result['stage3'] = {} for r in ('2', '3', '4'): res = sweep_operator(f'Lambda{r}', EXT[r], results) gate(f'Stage 3: M_{r} has no hyperexponential solution (0 candidates; every combination NOSOL by full column rank)', res['candidates'] == [], f"{res['ncombos_solved']} admissible combinations") gate(f'Stage 3: no completeness flag for M_{r}', res['flags'] == []) cs_sym, factors, rows, lams, lcomplete = indicial_tables(EXT[r]) ok_c, census = root_census(rows) gate(f'Stage 3: root census for M_{r} -- all indicial roots found at every factor; infinity complete', ok_c and lcomplete is True) result['stage3'][f'M_{r}'] = {'sweep': res, 'exponent_options': [{'factor': str(F.as_expr())[:40], 'deg': F.degree(), 'options': [str(e) for e in opts], 'complete': str(c)} for F, mult, v, k, opts, c, ind in rows], 'lam_options': [str(l) for l in lams], 'root_census': census} # ---------------- verdict ---------------- verdict = (nfail == 0) print(f'\nSUMMARY: {npass} PASS, {nfail} FAIL') if verdict: print('CERTIFICATE PASS (fail-closed): L5 at (2,3,5) has no right factor of order 1, 2, 3, or 4 over Q(xi)') else: print('CERTIFICATE FAIL: the lemma did NOT verify') result['checks'] = checks result['npass'], result['nfail'] = npass, nfail result['verdict'] = 'PASS' if verdict else 'FAIL' result['seconds'] = round(time.time() - t0, 1) json.dump(result, open(RESNAME, 'w'), indent=1) print(f'saved {RESNAME} ({result["seconds"]} s)') sys.stdout.flush() sys.stdout = sys.__stdout__ LOG.close() return 0 if verdict else 1 if __name__ == '__main__': try: sys.exit(main()) except SystemExit: raise except BaseException as e: # fail closed on any exception print(f'CERTIFICATE FAIL: exception {type(e).__name__}: {e}') sys.stdout.flush() sys.stdout = sys.__stdout__ sys.exit(2)