perf: integer-only forward elimination for gauss_solve (#72)

REFERENCES.md

@@ -40,14 +40,27 @@ f64 entry is decomposed into `mantissa × 2^exponent`, scaled to a common intege
 eliminated without any `BigRational` or GCD overhead.
 See `src/exact.rs` for the full architecture description.
 
-### f64 → BigRational conversion (`f64_to_bigrational`)
-
-`f64_to_bigrational` converts an f64 to an exact `BigRational` by decomposing the IEEE 754
-binary64 bit representation [9] into its sign, exponent, and significand fields.  Because
-every finite f64 is exactly `±m × 2^e` (where `m` is an integer), the rational can be
-constructed directly via `BigRational::new_raw` without GCD normalization — trailing zeros
-in the significand are stripped first so the fraction is already in lowest terms.  See
-Goldberg [10] for background on IEEE 754 representation and exact rational reconstruction.
+### Exact linear system solve (hybrid Bareiss / BigRational)
+
+`solve_exact()` and `solve_exact_f64()` share the BigInt core used for determinants.  Matrix
+and RHS entries are decomposed via IEEE 754 bit extraction [9] and scaled to a shared base
+`2^e_min` so the augmented system `(A | b)` becomes a `BigInt` matrix.  Forward elimination
+runs in `BigInt` using Bareiss fraction-free updates [7] — no `BigRational` and no GCD
+normalisation in the `O(D³)` phase.  The upper-triangular result is then lifted into
+`BigRational` for back-substitution, where fractions are inherent and the cost is only
+`O(D²)`.  Row swaps from first-non-zero pivoting are applied to both the matrix and the
+RHS; because power-of-two scaling is applied uniformly to both sides of `A x = b`, the
+solution is unchanged by the scale factor.
+
+### f64 → integer decomposition (`f64_decompose`)
+
+Both the determinant and solve paths convert f64 entries via `f64_decompose`, which extracts
+the IEEE 754 binary64 sign, unbiased exponent, and significand [9] and strips trailing zeros
+from the significand so `|x| = m · 2^e` with `m` odd.  The integer matrix is then assembled
+by shifting each mantissa left by `exp − e_min`, giving a GCD-free, Bareiss-ready starting
+point.  A one-shot wrapper `f64_to_bigrational` (used only in tests) packages the same
+decomposition into a single `BigRational`.  See Goldberg [10] for background on IEEE 754
+representation and exact rational reconstruction.
 
 ### LDL^T factorization (symmetric SPD/PSD)