1313//! `e - e_min`), and Bareiss elimination runs entirely in `BigInt`
1414//! arithmetic — no `BigRational`, no GCD, no denominator tracking.
1515//! The result is `(det_int, total_exp)` where `det = det_int × 2^(D × e_min)`.
16- //! `bareiss_det` wraps this with `bigint_exp_to_bigrational` to reconstruct
17- //! a reduced `BigRational`; `det_sign_exact` reads the sign directly from
18- //! `det_int` (the scale factor is always positive).
16+ //! `det_exact` wraps this with `bigint_exp_to_bigrational` to reconstruct a
17+ //! reduced `BigRational`; `det_exact_f64` converts the same pair directly to
18+ //! `f64`; and `det_sign_exact` reads the sign directly from `det_int` (the
19+ //! scale factor is always positive).
1920//!
2021//! `det_sign_exact` adds a two-stage adaptive-precision optimisation inspired
2122//! by Shewchuk's robust geometric predicates:
@@ -142,7 +143,8 @@ fn bigint_exp_to_bigrational(mut value: BigInt, mut exp: i32) -> BigRational {
142143 let exp_abs = exp. unsigned_abs ( ) ;
143144 let reduce = tz. min ( u64:: from ( exp_abs) ) ;
144145 value >>= reduce;
145- let reduce = u32:: try_from ( reduce) . unwrap_or ( u32:: MAX ) ;
146+ #[ allow( clippy:: cast_possible_truncation) ]
147+ let reduce = reduce as u32 ;
146148 let remaining_abs = exp_abs - reduce;
147149 exp = match remaining_abs {
148150 0 => 0 ,
@@ -158,6 +160,31 @@ fn bigint_exp_to_bigrational(mut value: BigInt, mut exp: i32) -> BigRational {
158160 }
159161}
160162
163+ /// Convert a `BigInt × 2^exp` determinant pair to finite `f64` without first
164+ /// reducing a public `BigRational` determinant value.
165+ fn bigint_exp_to_finite_f64 ( value : BigInt , exp : i32 ) -> Result < f64 , LaError > {
166+ if value == BigInt :: from ( 0 ) {
167+ return Ok ( 0.0 ) ;
168+ }
169+
170+ let exact = if exp >= 0 {
171+ BigRational :: new_raw ( value << exp. cast_unsigned ( ) , BigInt :: from ( 1u32 ) )
172+ } else {
173+ BigRational :: new_raw ( value, BigInt :: from ( 1u32 ) << exp. unsigned_abs ( ) )
174+ } ;
175+
176+ let Some ( val) = exact. to_f64 ( ) else {
177+ cold_path ( ) ;
178+ return Err ( LaError :: Overflow { index : None } ) ;
179+ } ;
180+ if val. is_finite ( ) {
181+ Ok ( val)
182+ } else {
183+ cold_path ( ) ;
184+ Err ( LaError :: Overflow { index : None } )
185+ }
186+ }
187+
161188// -----------------------------------------------------------------------
162189// Shared integer-Bareiss primitives
163190// -----------------------------------------------------------------------
@@ -396,11 +423,11 @@ fn bareiss_det_int_finite<const D: usize>(m: &FiniteMatrix<D>) -> Result<(BigInt
396423 // det(original) = det_int × 2^(D × e_min)
397424 let Ok ( d_i32) = i32:: try_from ( D ) else {
398425 cold_path ( ) ;
399- return Err ( LaError :: unsupported_dimension ( D , i32 :: MAX as usize ) ) ;
426+ return Err ( LaError :: determinant_scale_overflow ( D , e_min ) ) ;
400427 } ;
401428 let Some ( total_exp) = e_min. checked_mul ( d_i32) else {
402429 cold_path ( ) ;
403- return Err ( LaError :: Overflow { index : None } ) ;
430+ return Err ( LaError :: determinant_scale_overflow ( D , e_min ) ) ;
404431 } ;
405432
406433 Ok ( ( det_int, total_exp) )
@@ -489,21 +516,15 @@ impl<const D: usize> FiniteMatrix<D> {
489516 /// Exact determinant converted to a finite `f64`.
490517///
491518/// # Errors
519+ /// Returns [`LaError::DeterminantScaleOverflow`] if determinant scaling
520+ /// overflows the internal exponent representation.
521+ ///
492522/// Returns [`LaError::Overflow`] if the exact determinant cannot be
493523/// represented as a finite `f64`.
494524#[ inline]
495525 fn det_exact_f64 ( & self ) -> Result < f64 , LaError > {
496- let exact = self . det_exact ( ) ?;
497- let Some ( val) = exact. to_f64 ( ) else {
498- cold_path ( ) ;
499- return Err ( LaError :: Overflow { index : None } ) ;
500- } ;
501- if val. is_finite ( ) {
502- Ok ( val)
503- } else {
504- cold_path ( ) ;
505- Err ( LaError :: Overflow { index : None } )
506- }
526+ let ( det_int, total_exp) = bareiss_det_int_finite ( self ) ?;
527+ bigint_exp_to_finite_f64 ( det_int, total_exp)
507528 }
508529
509530 /// Exact linear solve for finite inputs.
@@ -548,6 +569,9 @@ impl<const D: usize> FiniteMatrix<D> {
548569/// Returns [`LaError::NonFinite`] if a direct determinant or error-bound
549570/// computation detects a non-finite condition that is not an inconclusive
550571/// scalar overflow.
572+ ///
573+ /// Returns [`LaError::DeterminantScaleOverflow`] if determinant scaling
574+ /// overflows the internal exponent representation.
551575#[ inline]
552576 fn det_sign_exact ( & self ) -> Result < i8 , LaError > {
553577 match ( self . det_direct ( ) , self . det_errbound ( ) ) {
@@ -607,11 +631,8 @@ impl<const D: usize> Matrix<D> {
607631/// # Errors
608632/// Returns [`LaError::NonFinite`] if stored matrix entries are NaN or infinity.
609633///
610- /// Returns [`LaError::Overflow`] if determinant scaling overflows the internal
611- /// exponent representation.
612- ///
613- /// Returns [`LaError::UnsupportedDimension`] if `D` cannot be represented in
614- /// the internal determinant exponent calculation.
634+ /// Returns [`LaError::DeterminantScaleOverflow`] if determinant scaling
635+ /// overflows the internal exponent representation.
615636#[ inline]
616637 pub fn det_exact ( & self ) -> Result < BigRational , LaError > {
617638 FiniteMatrix :: new ( * self ) ?. det_exact ( )
@@ -621,10 +642,12 @@ impl<const D: usize> Matrix<D> {
621642///
622643/// Requires the `exact` Cargo feature.
623644///
624- /// Computes the exact [`BigRational`] determinant via [`det_exact`](Self::det_exact)
625- /// and converts it to the nearest `f64`. This is useful when you want the
626- /// most accurate f64 determinant possible without committing to `BigRational`
627- /// in your downstream code.
645+ /// Computes the exact determinant with the same integer Bareiss core used by
646+ /// [`det_exact`](Self::det_exact), then converts the exact scaled integer
647+ /// result to the nearest `f64` without first materializing the public
648+ /// [`BigRational`] determinant. This is useful when you want the most accurate
649+ /// f64 determinant possible without committing to `BigRational` in your
650+ /// downstream code.
628651///
629652/// # Examples
630653/// ```
@@ -641,8 +664,10 @@ impl<const D: usize> Matrix<D> {
641664/// # Errors
642665/// Returns [`LaError::NonFinite`] if stored matrix entries are NaN or infinity.
643666///
644- /// Returns [`LaError::Overflow`] if determinant scaling overflows the internal
645- /// exponent representation or if the exact determinant is too large to
667+ /// Returns [`LaError::DeterminantScaleOverflow`] if determinant scaling
668+ /// overflows the internal exponent representation.
669+ ///
670+ /// Returns [`LaError::Overflow`] if the exact determinant is too large to
646671/// represent as a finite `f64`.
647672#[ inline]
648673 pub fn det_exact_f64 ( & self ) -> Result < f64 , LaError > {
@@ -778,7 +803,9 @@ impl<const D: usize> Matrix<D> {
778803///
779804/// # Errors
780805/// Returns [`LaError::NonFinite`] if stored matrix entries are NaN or infinity.
781- /// This exact sign path has no additional runtime errors for finite matrices.
806+ ///
807+ /// Returns [`LaError::DeterminantScaleOverflow`] if determinant scaling
808+ /// overflows the internal exponent representation.
782809#[ inline]
783810 pub fn det_sign_exact ( & self ) -> Result < i8 , LaError > {
784811 FiniteMatrix :: new ( * self ) ?. det_sign_exact ( )
0 commit comments