11//! Benchmarks for exact arithmetic operations.
22//!
33//! These benchmarks measure the performance of the `exact` feature's
4- //! arbitrary-precision methods across dimensions D=2..5 (the primary
5- //! target for geometric predicates).
4+ //! arbitrary-precision methods. They are organised into two classes:
5+ //!
6+ //! 1. **General-case benches** (`exact_d{2..5}`) — a single
7+ //! well-conditioned diagonally-dominant matrix per dimension. These
8+ //! measure typical-case performance and track regressions against a
9+ //! reproducible input.
10+ //! 2. **Adversarial / extreme-input benches** — matrices chosen to
11+ //! stress specific corners of the exact-arithmetic pipeline:
12+ //! near-singularity (forces the Bareiss fallback), large f64 entries
13+ //! (stresses intermediate `BigInt` growth), and Hilbert-style
14+ //! ill-conditioning (wide range of `(mantissa, exponent)` pairs in
15+ //! the `f64_decompose → BigInt` path). These measure tail behaviour
16+ //! that fixed well-conditioned inputs miss and provide stronger
17+ //! empirical evidence for `docs/PERFORMANCE.md`.
618
7- use criterion:: Criterion ;
19+ use criterion:: { BenchmarkGroup , Criterion , measurement :: WallTime } ;
820use la_stack:: { Matrix , Vector } ;
921use pastey:: paste;
1022use std:: hint:: black_box;
@@ -47,8 +59,13 @@ fn make_vector_array<const D: usize>() -> [f64; D] {
4759}
4860
4961/// Near-singular matrix: base singular matrix + tiny perturbation.
50- /// This forces the exact Bareiss fallback in `det_sign_exact` (the fast
51- /// f64 filter cannot resolve the sign).
62+ ///
63+ /// The base `[[1,2,3],[4,5,6],[7,8,9]]` is exactly singular; adding
64+ /// `2^-50` to entry (0,0) makes `det = -3 × 2^-50 ≠ 0`. The f64 filter
65+ /// in `det_sign_exact` cannot resolve this sign, so Bareiss is forced;
66+ /// `solve_exact` is the primary use case for near-degenerate inputs
67+ /// (exact circumcenter etc.) and exercises the largest intermediate
68+ /// `BigInt` values in the hybrid solve.
5269#[ inline]
5370fn near_singular_3x3 ( ) -> Matrix < 3 > {
5471 let perturbation = f64:: from_bits ( 0x3CD0_0000_0000_0000 ) ; // 2^-50
@@ -59,6 +76,97 @@ fn near_singular_3x3() -> Matrix<3> {
5976 ] )
6077}
6178
79+ /// Large-entry 3×3: strictly diagonally-dominant matrix with diagonal
80+ /// entries near `f64::MAX / 2` and ones elsewhere.
81+ ///
82+ /// Each big entry decomposes into a 53-bit mantissa with exponent `~970`;
83+ /// the unit off-diagonals have exponent `0`, so the shared `e_min = 0`
84+ /// shift in `component_to_bigint` produces `BigInt`s of `~1023` bits for
85+ /// the diagonal and small integers elsewhere. Bareiss fraction-free
86+ /// updates then multiply these together, stressing the big-integer
87+ /// multiply and allocator along the full `O(D³)` elimination phase. The
88+ /// matrix is non-singular (det ≈ `big³`) so both `det_*` and `solve_*`
89+ /// paths complete.
90+ #[ inline]
91+ fn large_entries_3x3 ( ) -> Matrix < 3 > {
92+ let big = f64:: MAX / 2.0 ;
93+ Matrix :: < 3 > :: from_rows ( [ [ big, 1.0 , 1.0 ] , [ 1.0 , big, 1.0 ] , [ 1.0 , 1.0 , big] ] )
94+ }
95+
96+ /// Hilbert matrix `H[i][j] = 1 / (i + j + 1)`.
97+ ///
98+ /// Most entries (`1/3`, `1/5`, `1/6`, `1/7`, …) are non-terminating in
99+ /// binary, so every cell has a distinct 53-bit mantissa and a small
100+ /// negative exponent. `f64_decompose` therefore produces a wide mix of
101+ /// `(mantissa, exponent)` pairs with no shared power-of-two factors,
102+ /// and the scaling shift to the common `e_min` yields `BigInt` values
103+ /// of varied bit-lengths — a different kind of adversarial input from
104+ /// the large-entries case. Hilbert matrices are also classically
105+ /// ill-conditioned (condition number grows exponentially with D), so
106+ /// they are a realistic stand-in for the near-degenerate geometric
107+ /// predicate inputs that motivate exact arithmetic.
108+ #[ inline]
109+ #[ allow( clippy:: cast_precision_loss) ]
110+ fn hilbert < const D : usize > ( ) -> Matrix < D > {
111+ let mut rows = [ [ 0.0 ; D ] ; D ] ;
112+ let mut r = 0 ;
113+ while r < D {
114+ let mut c = 0 ;
115+ while c < D {
116+ rows[ r] [ c] = 1.0 / ( ( r + c + 1 ) as f64 ) ;
117+ c += 1 ;
118+ }
119+ r += 1 ;
120+ }
121+ Matrix :: < D > :: from_rows ( rows)
122+ }
123+
124+ /// Populate a Criterion group with the four headline exact-arithmetic
125+ /// benches on a single `(matrix, rhs)` pair: `det_sign_exact`,
126+ /// `det_exact`, `solve_exact`, `solve_exact_f64`.
127+ ///
128+ /// Used by every adversarial-input group so each one measures the same
129+ /// operations, making the resulting tables directly comparable.
130+ fn bench_extreme_group < const D : usize > (
131+ group : & mut BenchmarkGroup < ' _ , WallTime > ,
132+ m : Matrix < D > ,
133+ rhs : Vector < D > ,
134+ ) {
135+ group. bench_function ( "det_sign_exact" , |bencher| {
136+ bencher. iter ( || {
137+ let sign = black_box ( m)
138+ . det_sign_exact ( )
139+ . expect ( "finite matrix entries" ) ;
140+ black_box ( sign) ;
141+ } ) ;
142+ } ) ;
143+
144+ group. bench_function ( "det_exact" , |bencher| {
145+ bencher. iter ( || {
146+ let det = black_box ( m) . det_exact ( ) . expect ( "finite matrix entries" ) ;
147+ black_box ( det) ;
148+ } ) ;
149+ } ) ;
150+
151+ group. bench_function ( "solve_exact" , |bencher| {
152+ bencher. iter ( || {
153+ let x = black_box ( m)
154+ . solve_exact ( black_box ( rhs) )
155+ . expect ( "non-singular matrix with finite entries" ) ;
156+ let _ = black_box ( x) ;
157+ } ) ;
158+ } ) ;
159+
160+ group. bench_function ( "solve_exact_f64" , |bencher| {
161+ bencher. iter ( || {
162+ let x = black_box ( m)
163+ . solve_exact_f64 ( black_box ( rhs) )
164+ . expect ( "solution representable in f64" ) ;
165+ let _ = black_box ( x) ;
166+ } ) ;
167+ } ) ;
168+ }
169+
62170macro_rules! gen_exact_benches_for_dim {
63171 ( $c: expr, $d: literal) => {
64172 paste! { {
@@ -72,7 +180,7 @@ macro_rules! gen_exact_benches_for_dim {
72180 bencher. iter( || {
73181 let det = black_box( a)
74182 . det( la_stack:: DEFAULT_PIVOT_TOL )
75- . expect( "should not fail " ) ;
183+ . expect( "diagonally dominant matrix is non-singular " ) ;
76184 black_box( det) ;
77185 } ) ;
78186 } ) ;
@@ -87,39 +195,47 @@ macro_rules! gen_exact_benches_for_dim {
87195 // === det_exact (BigRational result) ===
88196 [ <group_d $d>] . bench_function( "det_exact" , |bencher| {
89197 bencher. iter( || {
90- let det = black_box( a) . det_exact( ) . expect( "should not fail " ) ;
198+ let det = black_box( a) . det_exact( ) . expect( "finite matrix entries " ) ;
91199 black_box( det) ;
92200 } ) ;
93201 } ) ;
94202
95203 // === det_exact_f64 (exact → f64) ===
96204 [ <group_d $d>] . bench_function( "det_exact_f64" , |bencher| {
97205 bencher. iter( || {
98- let det = black_box( a) . det_exact_f64( ) . expect( "should not fail" ) ;
206+ let det = black_box( a)
207+ . det_exact_f64( )
208+ . expect( "det representable in f64" ) ;
99209 black_box( det) ;
100210 } ) ;
101211 } ) ;
102212
103213 // === det_sign_exact (adaptive: fast filter + exact fallback) ===
104214 [ <group_d $d>] . bench_function( "det_sign_exact" , |bencher| {
105215 bencher. iter( || {
106- let sign = black_box( a) . det_sign_exact( ) . expect( "should not fail" ) ;
216+ let sign = black_box( a)
217+ . det_sign_exact( )
218+ . expect( "finite matrix entries" ) ;
107219 black_box( sign) ;
108220 } ) ;
109221 } ) ;
110222
111223 // === solve_exact (BigRational result) ===
112224 [ <group_d $d>] . bench_function( "solve_exact" , |bencher| {
113225 bencher. iter( || {
114- let x = black_box( a) . solve_exact( black_box( rhs) ) . expect( "should not fail" ) ;
226+ let x = black_box( a)
227+ . solve_exact( black_box( rhs) )
228+ . expect( "diagonally dominant matrix is non-singular" ) ;
115229 black_box( x) ;
116230 } ) ;
117231 } ) ;
118232
119233 // === solve_exact_f64 (exact → f64) ===
120234 [ <group_d $d>] . bench_function( "solve_exact_f64" , |bencher| {
121235 bencher. iter( || {
122- let x = black_box( a) . solve_exact_f64( black_box( rhs) ) . expect( "should not fail" ) ;
236+ let x = black_box( a)
237+ . solve_exact_f64( black_box( rhs) )
238+ . expect( "solution representable in f64" ) ;
123239 black_box( x) ;
124240 } ) ;
125241 } ) ;
@@ -140,25 +256,50 @@ fn main() {
140256 gen_exact_benches_for_dim ! ( & mut c, 5 ) ;
141257 }
142258
143- // Near-singular 3×3: forces Bareiss fallback in det_sign_exact.
259+ // === Adversarial / extreme-input groups ===
260+ //
261+ // Each group runs the same four exact-arithmetic benches
262+ // (`det_sign_exact`, `det_exact`, `solve_exact`, `solve_exact_f64`)
263+ // via `bench_extreme_group`, so the resulting tables are directly
264+ // comparable across input classes.
265+
266+ // Near-singular 3×3: forces Bareiss fallback in det_sign_exact and
267+ // exercises the largest intermediate BigInt values in solve_exact
268+ // (the primary motivating use case for exact solve).
144269 {
145- let m = near_singular_3x3 ( ) ;
146270 let mut group = c. benchmark_group ( "exact_near_singular_3x3" ) ;
271+ bench_extreme_group (
272+ & mut group,
273+ near_singular_3x3 ( ) ,
274+ Vector :: < 3 > :: new ( [ 1.0 , 2.0 , 3.0 ] ) ,
275+ ) ;
276+ group. finish ( ) ;
277+ }
147278
148- group. bench_function ( "det_sign_exact" , |bencher| {
149- bencher. iter ( || {
150- let sign = black_box ( m) . det_sign_exact ( ) . expect ( "should not fail" ) ;
151- black_box ( sign) ;
152- } ) ;
153- } ) ;
279+ // Large-entry 3×3: diagonal entries near `f64::MAX / 2` stress
280+ // BigInt growth during Bareiss forward elimination.
281+ {
282+ let mut group = c. benchmark_group ( "exact_large_entries_3x3" ) ;
283+ bench_extreme_group (
284+ & mut group,
285+ large_entries_3x3 ( ) ,
286+ Vector :: < 3 > :: new ( [ 1.0 , 1.0 , 1.0 ] ) ,
287+ ) ;
288+ group. finish ( ) ;
289+ }
154290
155- group. bench_function ( "det_exact" , |bencher| {
156- bencher. iter ( || {
157- let det = black_box ( m) . det_exact ( ) . expect ( "should not fail" ) ;
158- black_box ( det) ;
159- } ) ;
160- } ) ;
291+ // Hilbert 4×4 and 5×5: classically ill-conditioned matrices whose
292+ // entries span many orders of magnitude in `(mantissa, exponent)`
293+ // space, exercising the f64 → BigInt scaling path.
294+ {
295+ let mut group = c. benchmark_group ( "exact_hilbert_4x4" ) ;
296+ bench_extreme_group ( & mut group, hilbert :: < 4 > ( ) , Vector :: < 4 > :: new ( [ 1.0 ; 4 ] ) ) ;
297+ group. finish ( ) ;
298+ }
161299
300+ {
301+ let mut group = c. benchmark_group ( "exact_hilbert_5x5" ) ;
302+ bench_extreme_group ( & mut group, hilbert :: < 5 > ( ) , Vector :: < 5 > :: new ( [ 1.0 ; 5 ] ) ) ;
162303 group. finish ( ) ;
163304 }
164305
0 commit comments