2222//!
2323//! ## Linear system solve
2424//!
25- //! `solve_exact` and `solve_exact_f64` solve `A x = b` using standard Gaussian
26- //! elimination with partial pivoting in `BigRational` arithmetic. Every finite
27- //! `f64` is exactly representable as a rational, so the result is provably
28- //! correct.
25+ //! `solve_exact` and `solve_exact_f64` solve `A x = b` using Gaussian
26+ //! elimination with first-non-zero pivoting in `BigRational` arithmetic.
27+ //! Since all arithmetic is exact, any non-zero pivot gives the correct result
28+ //! (there is no numerical stability concern). Every finite `f64` is exactly
29+ //! representable as a rational, so the result is provably correct.
2930
3031use num_bigint:: { BigInt , Sign } ;
3132use num_rational:: BigRational ;
@@ -101,7 +102,7 @@ fn bareiss_det<const D: usize>(m: &Matrix<D>) -> BigRational {
101102 let mut sign: i8 = 1 ;
102103
103104 for k in 0 ..D {
104- // Partial pivoting: find non-zero entry in column k at or below row k.
105+ // Find first non-zero entry in column k at or below row k.
105106 if a[ k] [ k] == zero {
106107 let mut found = false ;
107108 for i in ( k + 1 ) ..D {
@@ -136,8 +137,12 @@ fn bareiss_det<const D: usize>(m: &Matrix<D>) -> BigRational {
136137 if sign < 0 { -det } else { det. clone ( ) }
137138}
138139
139- /// Solve `A x = b` using Gaussian elimination with partial pivoting in
140- /// `BigRational` arithmetic.
140+ /// Solve `A x = b` using Gaussian elimination with first-non-zero pivoting
141+ /// in `BigRational` arithmetic.
142+ ///
143+ /// Since all arithmetic is exact, any non-zero pivot gives the correct result
144+ /// (no numerical stability concern). This matches the pivoting strategy used
145+ /// by `bareiss_det`.
141146///
142147/// Returns the exact solution as a `Vec<BigRational>` of length `D`.
143148/// Returns `Err(LaError::Singular)` if the matrix is exactly singular.
@@ -159,7 +164,7 @@ fn gauss_solve<const D: usize>(m: &Matrix<D>, b: &Vector<D>) -> Result<Vec<BigRa
159164 aug. push ( row) ;
160165 }
161166
162- // Forward elimination with partial pivoting.
167+ // Forward elimination with first-non-zero pivoting.
163168 for k in 0 ..D {
164169 // Find first non-zero pivot in column k at or below row k.
165170 if aug[ k] [ k] == zero {
@@ -189,8 +194,7 @@ fn gauss_solve<const D: usize>(m: &Matrix<D>, b: &Vector<D>) -> Result<Vec<BigRa
189194
190195 // Back-substitution.
191196 let mut x: Vec < BigRational > = vec ! [ zero; D ] ;
192- for ii in 0 ..D {
193- let i = D - 1 - ii;
197+ for i in ( 0 ..D ) . rev ( ) {
194198 let mut sum = aug[ i] [ D ] . clone ( ) ;
195199 for j in ( i + 1 ) ..D {
196200 sum -= & aug[ i] [ j] * & x[ j] ;
0 commit comments