perf: integer-only forward elimination for gauss_solve (BigInt hybrid)

[acgetchell]

Summary

Apply the integer-only Bareiss approach from #64 to the forward elimination phase of gauss_solve, while keeping rational back-substitution. This hybrid avoids BigRational GCD overhead during the most expensive phase of the solve.

Current State

gauss_solve operates entirely in BigRational: D² matrix entries + D RHS entries are converted via f64_to_bigrational, then standard Gaussian elimination with rational division + back-substitution produces the exact solution.

With #64 merged, the determinant path runs ~16–40x faster using pure BigInt. The solve path is unchanged (solve_exact benchmarks show 1.00x).

Proposed Changes

Hybrid integer/rational solve

1. Forward elimination in BigInt: decompose matrix and RHS into mantissa × 2^exponent, find e_min, scale to integers, run elimination with exact integer division (same approach as bareiss_det_int)
2. Back-substitution in BigRational: after forward elimination, convert the upper-triangular BigInt matrix + RHS back to BigRational for back-substitution (which inherently produces fractions)

Design considerations

• The RHS vector has its own exponents — need to track a separate e_min_rhs or a unified minimum across matrix + RHS
• Row swaps must keep the RHS in sync (same as current code)
• The pivot division in forward elimination is factor = mat[i][k] / pivot — in the integer path this becomes the Bareiss-style (pivot * a[i][j] - a[i][k] * a[k][j]) / prev_pivot pattern
• Back-substitution: x[i] = (rhs[i] - sum) / mat[i][i] — this division is generally not exact in integers, so BigRational is needed here

Benefits

• Eliminates GCD normalization during forward elimination (the O(D³) phase)
• Back-substitution is only O(D²) — rational overhead is acceptable there
• Expected speedup: significant for D=4–5 where forward elimination dominates

Implementation Notes

• Depends on: perf: integer-only Bareiss determinant (BigInt path) #64 (integer Bareiss, merged)
• Reuse f64_decompose and the exponent-scaling pattern from bareiss_det_int
• The current gauss_solve uses standard Gaussian elimination (not Bareiss) — the hybrid would switch forward elimination to Bareiss-style fraction-free form
• All existing solve_exact and solve_exact_f64 tests must continue to pass
• Add benchmarks comparing before/after for D=2–5

[acgetchell]

Blocks #80 (benchmark coverage should measure the improved solve path).