diff --git a/formal/SilverSight/PIST/CharPoly.lean b/formal/SilverSight/PIST/CharPoly.lean index 7acf58b6..2e3b1f9f 100644 --- a/formal/SilverSight/PIST/CharPoly.lean +++ b/formal/SilverSight/PIST/CharPoly.lean @@ -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) let result := (List.range n).foldl (fun acc i => let c := coeffs.getD i 0 - add (mul acc x) (ofRawInt c)) x + add (mul acc x) (ofInt c)) one result 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 c := coeffs.getD i 0 - let coeff := ofRawInt ((n - i : Int) * c) - add (mul acc x) coeff) (ofRawInt (n : Int)) + let coeff := ofInt ((n - 1 - i : Int) * c) + add (mul acc x) coeff) (ofInt (n : Int)) result -- Initial guess: use the Frobenius norm as a rough bound 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. #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 end SilverSight.PIST.CharPoly