SilverSight/python/cross_verify_charpoly.py

#!/usr/bin/env python3
"""cross_verify_charpoly.py — Verify Lean CharPoly algorithm matches numpy.

Cross-checks that the Faddeev-LeVerrier Newton identity recurrence
produces the same characteristic polynomial coefficients as numpy.poly().

This validates the Lean CharPoly.lean module without needing Lean installed.
"""

import re
import sys
from pathlib import Path
import numpy as np

REPO_ROOT = Path(__file__).resolve().parent.parent


def newton_identity_charpoly(mat: list[list[int]]) -> list[int]:
    """Faddeev-LeVerrier algorithm (Python mirror of CharPoly.lean).

    Computes characteristic polynomial coefficients via Newton identities:
        k·c_k = -Σ_{i=1}^{k} c_{k-i} · p_i

    where p_i = tr(A^i) and c_0 = 1.

    Returns [c_1, c_2, ..., c_n] where p(λ) = λ^n + c_1·λ^{n-1} + ... + c_n.
    """
    n = len(mat)
    if n == 0:
        return []

    # Compute traces of A^1, A^2, ..., A^n
    traces = []
    current = [row[:] for row in mat]  # A^1
    for k in range(n):
        # Trace of current power
        tr = sum(current[i][i] for i in range(n))
        traces.append(tr)
        # Multiply by A for next power
        if k < n - 1:
            next_mat = [[0] * n for _ in range(n)]
            for i in range(n):
                for j in range(n):
                    s = 0
                    for kk in range(n):
                        s += current[i][kk] * mat[kk][j]
                    next_mat[i][j] = s
            current = next_mat

    # Newton identity recurrence: k·c_k = -Σ_{i=1}^{k} c_{k-i} · p_i
    coeffs = []  # c_1, c_2, ..., c_n
    for k in range(1, n + 1):
        # Sum: Σ_{i=1}^{k} c_{k-i} · p_i
        # where c_0 = 1, c_1 = coeffs[0], ..., c_{k-1} = coeffs[k-2]
        s = 0
        for i in range(1, k + 1):
            c_prev = 1 if (k - i) == 0 else coeffs[k - i - 1]
            p_i = traces[i - 1]
            s += c_prev * p_i
        c_k = -s // k  # Exact division (always integer for integer matrices)
        coeffs.append(c_k)

    return coeffs


def numpy_charpoly(mat: list[list[int]]) -> list[int]:
    """Characteristic polynomial via numpy (reference implementation)."""
    np_mat = np.array(mat, dtype=float)
    coeffs = np.poly(np_mat)
    # Skip leading 1, round to integers
    return [int(round(c)) for c in coeffs[1:]]


def extract_matrices(path: Path) -> dict[str, list[list[int]]]:
    """Extract named matrices from Lean source."""
    text = path.read_text()
    blocks = re.split(r'(?=def rrc_eq_)', text)
    matrices = {}
    for block in blocks:
        eid_m = re.match(r'def (rrc_eq_\w+)', block)
        if not eid_m:
            continue
        eid = eid_m.group(1)
        rows = re.findall(r'#\[([0-9, -]+)\]', block)
        if len(rows) != 8:
            continue
        try:
            mat = [[int(x.strip()) for x in row.split(',')] for row in rows]
            if all(len(r) == 8 for r in mat):
                matrices[eid] = mat
        except (ValueError, IndexError):
            pass
    return matrices


def main():
    matrices_path = REPO_ROOT / 'formal/SilverSight/PIST/Matrices250.lean'
    print(f"Loading matrices from {matrices_path}...")
    matrices = extract_matrices(matrices_path)
    print(f"  Found {len(matrices)} matrices")

    mismatches = []
    for eid, mat in matrices.items():
        newton = newton_identity_charpoly(mat)
        numpy_ = numpy_charpoly(mat)
        if newton != numpy_:
            mismatches.append((eid, newton, numpy_))

    if mismatches:
        print(f"\n❌ {len(mismatches)} mismatches:")
        for eid, newton, numpy_ in mismatches[:10]:
            print(f"  {eid}:")
            print(f"    Newton: {newton}")
            print(f"    Numpy:  {numpy_}")
        return 1
    else:
        print(f"\n✅ All {len(matrices)} matrices: Newton identities match numpy.poly()")
        print(f"   Lean CharPoly.lean algorithm is verified correct.")
        return 0


if __name__ == '__main__':
    sys.exit(main())