#!/usr/bin/env python3 """Generate kernel-checked fraction-free integer moment certificates. The Python calculation is untrusted: generated Lean proofs check every Schur identity and every positive pivot using ordinary kernel reduction. """ from __future__ import annotations import argparse from hashlib import sha256 from collections import defaultdict from pathlib import Path import check_certificate as certificate ROOT = Path(__file__).resolve().parents[1] CHECK = False CERTIFICATE_SHA256 = '2417f05c4651a380467d1c6ca7d10efcf7119790d5d3c86af939e35a5843779b' def emit(path, content): if CHECK: if not path.exists() or path.read_text() != content: raise ValueError(f'generated Lean source differs: {path.name}') else: path.write_text(content) HEADER = '''module public import PersistentSOS.PSD @[expose] public section /-! # Exact moment block certificate Generated by `scripts/generate_moments.py`; all identities below are checked by the Lean kernel with integer arithmetic. No external verifier is trusted. -/ namespace PersistentSOS set_option maxRecDepth 100000 set_option maxHeartbeats 0 ''' def vector(values): return '![' + ', '.join(map(str, values)) + ']' def matrix(values): return '![' + ',\n '.join(vector(row) for row in values) + ']' def schur_stages(a): stages = [a] previous = 1 while a: pivot = a[0][0] if pivot == 0: certificate.require(all(v == 0 for row in a for v in row), 'zero pivot in nonzero matrix') break certificate.require(pivot > 0, 'nonpositive Schur pivot') b = [] for i in range(1, len(a)): row = [] for j in range(1, len(a)): value, remainder = divmod(a[i][j] * pivot - a[i][0] * a[0][j], previous) certificate.require(remainder == 0, 'nonexact Schur division') row.append(value) b.append(row) stages.append(b) previous, a = pivot, b return stages def generate_block(ch, basis, moment): original = [[moment(tuple(sorted(u + v))) for v in basis] for u in basis] stages = schur_stages(original) out = [HEADER] terminal_zero = bool(stages[-1]) for k, a in enumerate(stages if terminal_zero else stages[:-1]): n = len(a) out.append(f'def momentBlock{ch}Stage{k} : Fin {n} → Fin {n} → ℤ :=\n {matrix(a)}\n\n') if terminal_zero: k = len(stages)-1 out.append(f'theorem momentBlock{ch}Stage{k}_psd : IsPSD (fun i j => (momentBlock{ch}Stage{k} i j : ℝ)) := by\n intro x\n simp [quadraticForm, momentBlock{ch}Stage{k}]\n\n') else: out.append(f'def momentBlock{ch}Stage{len(original)} : Fin 0 → Fin 0 → ℤ := fun i => Fin.elim0 i\n\n') for k in reversed(range(len(stages)-1)): a = stages[k] n = len(a) name = f'momentBlock{ch}Stage{k}' tail = f'momentBlock{ch}Stage{k+1}' previous = stages[k-1][0][0] if k else 1 out.append(f'theorem {name}_psd : IsPSD (fun i j => ({name} i j : ℝ)) := by\n') out.append(f' apply integer_schur_step {name} {tail} ({previous})\n') out.append(' · decide\n · decide\n') out.append(' · decide +kernel\n · decide +kernel\n') if n == 1: out.append(' · exact quadraticForm_zero _\n\n') else: out.append(f' · exact {tail}_psd\n\n') out.append('end PersistentSOS\n') path = ROOT / 'PersistentSOS' / f'MomentBlock{ch}.lean' emit(path, ''.join(out)) print(path.name, len(original), path.stat().st_size) def main(): moment = certificate.functional(certificate.load_certificate()) blocks = defaultdict(list) for degree in range(3): for monomial in certificate.monomials(degree): if monomial != (0,0): blocks[certificate.character(monomial)].append(monomial) for ch, basis in blocks.items(): if ch != 0: generate_block(ch, basis, moment) generate_block("Full0", blocks[0] + [(0,0)], moment) def generate_data(): moment = certificate.functional(certificate.load_certificate()) blocks = defaultdict(list) for degree in range(3): for monomial in certificate.monomials(degree): if monomial != (0,0): blocks[certificate.character(monomial)].append(monomial) blocks[0].append((0,0)) data = {} for degree in range(5): for monomial in certificate.monomials(degree): value = moment(monomial) if value: data[sum(5**i for i in monomial)] = value def tree(rows, indent): if not rows: return '0' if len(rows) == 1: k,v = rows[0] return f'if k = {k} then ({v}) else 0' mid = len(rows)//2 pivot = rows[mid][0] return f'if k < {pivot} then\n'+indent+' '+tree(rows[:mid],indent+' ')+ '\n'+indent+'else\n'+indent+' '+tree(rows[mid:],indent+' ') out = ['module\n\npublic import PersistentSOS.MomentGeneral\n'] for c in range(8): out.append(f'public import PersistentSOS.MomentBlock{"Full0" if c == 0 else c}\n') out.append(''' /-! # Explicit rational moment data Every monomial occurrence uses one base-five exponent key. All degree-four moments are derived from the retained rational certificate. The finite Gram identities and coverage are proved below by kernel reduction. -/ @[expose] public section namespace PersistentSOS set_option maxRecDepth 100000 set_option maxHeartbeats 0 /-- Integer numerator table; unspecified keys have value zero. -/ def momentNumerator (k : ℕ) : ℤ := ''') out.append(tree(sorted(data.items()),' ')+'\n\n') out.append('def blockSize : Fin 8 → ℕ\n') for c in range(8): out.append(f' | {c} => {len(blocks[c])}\n') out.append('\ndef basisPair : (c : Fin 8) → Fin (blockSize c) → Fin 15 × Fin 15\n') for c in range(8): pairs = [(tuple([0]*(2-len(m)))+tuple(i+1 for i in m)) for m in blocks[c]] out.append(f' | {c} => {vector(pairs)}\n') out.append('\ndef blockNumerator : (c : Fin 8) → Fin (blockSize c) → Fin (blockSize c) → ℤ\n') for c in range(8): out.append(f' | {c} => momentBlock{"Full0" if c==0 else c}Stage0\n') out.append('\ntheorem blockNumerator_psd (c : Fin 8) : IsPSD (fun i j => (blockNumerator c i j : ℝ)) := by\n fin_cases c\n') for c in range(8): out.append(f' · exact momentBlock{"Full0" if c==0 else c}Stage0_psd\n') out.append(''' /-- A zero index pads a monomial of degree less than two. -/ noncomputable def variableExponent (i : Fin 15) : Exponent := if h : i = 0 then 0 else Finsupp.single ⟨i.val - 1, by have hi := i.isLt have hn : i.val ≠ 0 := by simpa using h omega⟩ 1 def variableKey (i : Fin 15) : ℕ := if i = 0 then 0 else 5 ^ (i.val - 1) noncomputable def momentBasis (b : (c : Fin 8) × Fin (blockSize c)) : Exponent := variableExponent (basisPair b.1 b.2).1 + variableExponent (basisPair b.1 b.2).2 def basisKey (c : Fin 8) (i : Fin (blockSize c)) : ℕ := variableKey (basisPair c i).1 + variableKey (basisPair c i).2 @[simp] theorem exponentKey_variableExponent (i : Fin 15) : exponentKey (variableExponent i) = variableKey i := by simp only [variableExponent, variableKey] split_ifs <;> simp @[simp] theorem exponentKey_momentBasis (b : (c : Fin 8) × Fin (blockSize c)) : exponentKey (momentBasis b) = basisKey b.1 b.2 := by simp [momentBasis, basisKey] ''') for c in range(8): out.append(f'theorem moment_block_{c} : ∀ i j : Fin (blockSize {c}),\n momentNumerator (basisKey {c} i + basisKey {c} j) = blockNumerator {c} i j := by\n decide +kernel\n\n') out.append('theorem moment_block : ∀ c (i j : Fin (blockSize c)),\n momentNumerator (basisKey c i + basisKey c j) = blockNumerator c i j := by\n intro c\n fin_cases c\n') for c in range(8): out.append(f' · exact moment_block_{c}\n') out.append('\n') for c in range(8): out.append(f'theorem moment_cross_{c} : ∀ d : Fin 8, {c} ≠ d → ∀ (i : Fin (blockSize {c})) (j : Fin (blockSize d)),\n momentNumerator (basisKey {c} i + basisKey d j) = 0 := by\n decide +kernel\n\n') out.append('theorem moment_cross : ∀ c d : Fin 8, c ≠ d → ∀ (i : Fin (blockSize c)) (j : Fin (blockSize d)),\n momentNumerator (basisKey c i + basisKey d j) = 0 := by\n intro c\n fin_cases c\n') for c in range(8): out.append(f' · exact moment_cross_{c}\n') out.append(''' theorem moment_quadratic_numerator_sum : (∑ i : Fin 14, momentNumerator (2 * 5 ^ i.val)) = 1000000000 := by decide +kernel theorem moment_quartic_numerator_sum : (∑ i : Fin 14, ∑ j : Fin 14, momentNumerator (2 * 5 ^ i.val + 2 * 5 ^ j.val)) = 1000000000 := by decide +kernel end PersistentSOS ''') path = ROOT/'PersistentSOS'/'MomentData.lean' emit(path, ''.join(out)) print(path.name, len(data), path.stat().st_size) if __name__ == '__main__': parser = argparse.ArgumentParser(description=__doc__) parser.add_argument('--check', action='store_true', help='compare generated files without modifying them') CHECK = parser.parse_args().check certificate.require(sha256(certificate.CERTIFICATE.read_bytes()).hexdigest() == CERTIFICATE_SHA256, 'certificate file hash changed') free = certificate.load_certificate() moment = certificate.functional(free) certificate.check_sphere(moment) certificate.moment_matrices(moment) certificate.require(certificate.objective(moment) == certificate.VALUE, 'objective changed') main() generate_data()