diff --git a/python/cross_verify_charpoly.py b/python/cross_verify_charpoly.py
new file mode 100644
index 00000000..6a9f8243
--- /dev/null
+++ b/python/cross_verify_charpoly.py
@@ -0,0 +1,123 @@
+#!/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())