allaun/SilverSight
1
0
Fork 0
mirror of https://github.com/allaunthefox/SilverSight.git synced 2026-07-31 01:25:21 +00:00
Code Issues Projects Releases Packages Wiki Activity Actions

fix: CharPoly Newton solver — correct Horner evaluation

Browse source 
Two bugs: (1) evalPoly started with x instead of one (leading coeff),
giving degree n+1 polynomial. (2) evalDeriv used (n-i)*ci instead of
(n-1-i)*ci, off by one. Both coefficients used ofRawInt (wrong scale)
instead of ofInt (correct Q16_16 scale).

Result: spectral radius of 2x2 identity now 65408 (~99.8% of 65536).
Remaining error is Q16_16 underflow for double root (lambda-1)^2.

Build: 3309 jobs, 0 errors
This commit is contained in:
allaun 2026-07-05 08:11:31 -05:00
parent ddb0145f54
commit 04546db009
1 changed files with 8 additions and 5 deletions
Show stats Download patch file Download diff file Expand all files Collapse all files

13
formal/SilverSight/PIST/CharPoly.lean
View file

@ -113,14 +113,14 @@ def newtonLargestRoot (coeffs : Array Int) (maxIter : Nat := 50) : Q16_16 :=
-- Horner's method: p(x) = (...((x + c1)x + c2)x + ... + cn) -- Horner's method: p(x) = (...((x + c1)x + c2)x + ... + cn)
let result := (List.range n).foldl (fun acc i => let result := (List.range n).foldl (fun acc i =>
let c := coeffs.getD i 0 let c := coeffs.getD i 0
add (mul acc x) (ofRawInt c)) x add (mul acc x) (ofInt c)) one
result result
let evalDeriv (x : Q16_16) : Q16_16 := let evalDeriv (x : Q16_16) : Q16_16 :=
-- p'(x) = n*xn?1 + (n-1)*c1*xn?2 + ... + cn?1 -- p'(x) = n*x^{n-1} + (n-1)*c1*x^{n-2} + ... + 1*c_{n-1}
let result := (List.range (n - 1)).foldl (fun acc i => let result := (List.range (n - 1)).foldl (fun acc i =>
let c := coeffs.getD i 0 let c := coeffs.getD i 0
let coeff := ofRawInt ((n - i : Int) * c) let coeff := ofInt ((n - 1 - i : Int) * c)
add (mul acc x) coeff) (ofRawInt (n : Int)) add (mul acc x) coeff) (ofInt (n : Int))
result result
-- Initial guess: use the Frobenius norm as a rough bound -- Initial guess: use the Frobenius norm as a rough bound
let maxCoeff := coeffs.foldl (fun acc c => max acc (Int.natAbs c)) 0 let maxCoeff := coeffs.foldl (fun acc c => max acc (Int.natAbs c)) 0
@ -187,7 +187,10 @@ def exactSpectralRadius (n : Nat) (mat : Array (Array Int)) : Q16_16 :=
-- Test: trace of 2x2 identity is 2. -- Test: trace of 2x2 identity is 2.
#eval matrixTrace 2 #[#[1, 0], #[0, 1]] #eval matrixTrace 2 #[#[1, 0], #[0, 1]]
-- Test: exact spectral radius of 2x2 identity is 1.0 (65536 in Q16_16). -- Test: exact spectral radius of 2x2 identity is ~1.0.
-- Q16_16 Newton converges to 65408 (99.8% of 65536 = 1.0).
-- The remaining error is from Q16_16 underflow: (x-1)^2 underflows to 0
-- when |x-1| < 256/65536 (double root convergence limit).
#eval (exactSpectralRadius 2 #[#[1, 0], #[0, 1]]).toInt #eval (exactSpectralRadius 2 #[#[1, 0], #[0, 1]]).toInt
end SilverSight.PIST.CharPoly end SilverSight.PIST.CharPoly