la_stack/rational.rs
1#![forbid(unsafe_code)]
2
3//! Exact-input fixed-size matrices and vectors.
4//!
5//! [`RationalMatrix`] and [`RationalVector`] preserve caller-supplied
6//! [`BigRational`] coefficients without a binary64 round trip. Triangular solves
7//! (including row permutations) use direct exact substitution. Determinants
8//! and general solves clear denominators with a positive scale per row, then
9//! reuse the
10//! crate's fraction-free [`BigInt`] Bareiss backend. The positive row scales
11//! preserve determinant sign; determinant values divide by their product; and
12//! general solves apply the same row scale to the matrix and right-hand side.
13//! See `REFERENCES.md` \[7\] for Bareiss elimination and \[12\] for substitution and determinant
14//! multilinearity. The
15//! [row-clearing construction](https://github.com/acgetchell/la-stack/blob/main/docs/mathematical_basis.md#exact-arithmetic-over-rational-inputs)
16//! explains how these identities apply to canonical rational inputs.
17
18use std::array::from_fn;
19
20use num_bigint::{BigInt, Sign};
21use num_rational::BigRational;
22use num_traits::{One, Zero};
23
24use crate::exact::{det_big_int, solve_big_int};
25use crate::{DeterminantSign, LaError};
26
27/// Exact rational square matrix with compile-time dimension `D`.
28///
29/// Construction validates that every denominator is non-zero and canonicalizes
30/// every entry to lowest terms with a positive denominator. The private storage
31/// then carries those invariants, so determinant and solve methods do not repeat
32/// input validation. Unlike [`crate::Matrix`], entries are already exact rational
33/// values rather than finite binary64 values interpreted exactly.
34///
35/// Direct field construction is intentionally unavailable:
36///
37/// ```compile_fail
38/// use la_stack::{BigRational, RationalMatrix};
39///
40/// let _ = RationalMatrix::<1> {
41/// rows: [[BigRational::from_integer(1.into())]],
42/// };
43/// ```
44#[derive(Clone, Debug, Eq, PartialEq)]
45#[must_use]
46pub struct RationalMatrix<const D: usize> {
47 rows: [[BigRational; D]; D],
48}
49
50/// Exact rational vector with compile-time dimension `D`.
51///
52/// Construction validates that every denominator is non-zero and canonicalizes
53/// every entry to lowest terms with a positive denominator. Solutions returned
54/// by [`RationalMatrix::solve`] and [`crate::Matrix::solve_exact`] also use this
55/// type, making any later conversion to [`Vector`](crate::Vector) explicit through
56/// [`ExactF64Conversion`](crate::ExactF64Conversion).
57/// Use [`try_to_f64`](crate::ExactF64Conversion::try_to_f64) when every component
58/// must remain exact, or [`to_rounded_f64`](crate::ExactF64Conversion::to_rounded_f64)
59/// to opt into round-to-nearest, ties-to-even. Both conversions reject results
60/// that cannot be rounded to finite binary64 values.
61#[derive(Clone, Debug, Eq, PartialEq)]
62#[must_use]
63pub struct RationalVector<const D: usize> {
64 data: [BigRational; D],
65}
66
67impl<const D: usize> RationalMatrix<D> {
68 /// Try to create an exact matrix from row-major rational storage.
69 ///
70 /// Raw, non-reduced [`BigRational::new_raw`] values and negative
71 /// denominators are accepted, interpreted as their mathematical quotient,
72 /// and stored canonically. A raw zero denominator is not a rational value
73 /// and is rejected at this construction boundary.
74 ///
75 /// # Examples
76 /// ```
77 /// use la_stack::prelude::*;
78 ///
79 /// # fn main() -> Result<(), LaError> {
80 /// let matrix = RationalMatrix::<2>::try_from_rows([
81 /// [
82 /// BigRational::new(1.into(), 3.into()),
83 /// BigRational::from_integer(2.into()),
84 /// ],
85 /// [
86 /// BigRational::from_integer(1.into()),
87 /// BigRational::new(5.into(), 2.into()),
88 /// ],
89 /// ])?;
90 /// assert_eq!(
91 /// matrix.det(),
92 /// BigRational::new((-7).into(), 6.into())
93 /// );
94 /// # Ok(())
95 /// # }
96 /// ```
97 ///
98 /// # Errors
99 /// Returns [`LaError::NonFinite`] at the first matrix cell whose raw
100 /// rational denominator is zero.
101 pub fn try_from_rows(rows: [[BigRational; D]; D]) -> Result<Self, LaError> {
102 for (row_index, row) in rows.iter().enumerate() {
103 for (col_index, value) in row.iter().enumerate() {
104 if value.denom().sign() == Sign::NoSign {
105 return Err(LaError::non_finite_input_matrix(row_index, col_index));
106 }
107 }
108 }
109 Ok(Self {
110 rows: rows.map(|row| row.map(canonicalize_rational)),
111 })
112 }
113
114 /// Try to create an exact matrix by evaluating a function at every cell.
115 ///
116 /// The function is evaluated once per cell in row-major order.
117 ///
118 /// # Examples
119 /// ```
120 /// use la_stack::prelude::*;
121 ///
122 /// # fn main() -> Result<(), LaError> {
123 /// let diagonal = RationalMatrix::<3>::try_from_fn(|row, col| {
124 /// BigRational::from_integer(u8::from(row == col).into())
125 /// })?;
126 /// assert_eq!(diagonal.det_sign(), DeterminantSign::Positive);
127 /// # Ok(())
128 /// # }
129 /// ```
130 ///
131 /// # Errors
132 /// Returns [`LaError::NonFinite`] at the first generated cell whose raw
133 /// rational denominator is zero.
134 pub fn try_from_fn(
135 mut make_entry: impl FnMut(usize, usize) -> BigRational,
136 ) -> Result<Self, LaError> {
137 let rows = from_fn(|row| from_fn(|col| make_entry(row, col)));
138 Self::try_from_rows(rows)
139 }
140
141 /// Return the all-zero exact matrix.
142 pub fn zero() -> Self {
143 Self {
144 rows: from_fn(|_| from_fn(|_| BigRational::from_integer(BigInt::from(0)))),
145 }
146 }
147
148 /// Borrow the row-major exact storage.
149 #[must_use]
150 pub const fn as_rows(&self) -> &[[BigRational; D]; D] {
151 &self.rows
152 }
153
154 /// Consume the matrix and return its row-major exact storage.
155 #[must_use]
156 pub fn into_rows(self) -> [[BigRational; D]; D] {
157 self.rows
158 }
159
160 /// Borrow one entry, returning `None` for an out-of-bounds index.
161 #[must_use]
162 pub fn get(&self, row: usize, col: usize) -> Option<&BigRational> {
163 self.rows.get(row)?.get(col)
164 }
165
166 /// Replace one exact entry while preserving the canonical non-zero-
167 /// denominator invariant.
168 ///
169 /// As with [`try_from_rows`](Self::try_from_rows), non-reduced values and
170 /// negative denominators are accepted and canonicalized. Rejected indices
171 /// or denominators leave the matrix unchanged.
172 ///
173 /// # Examples
174 /// ```
175 /// use core::assert_matches;
176 /// use la_stack::prelude::*;
177 ///
178 /// # fn main() -> Result<(), LaError> {
179 /// let mut matrix = RationalMatrix::<2>::zero();
180 /// matrix.set(0, 1, BigRational::new_raw((-2).into(), (-4).into()))?;
181 /// let half = BigRational::new(1.into(), 2.into());
182 /// assert_eq!(matrix.get(0, 1), Some(&half));
183 ///
184 /// let before = matrix.clone();
185 /// assert_matches!(
186 /// matrix.set(0, 1, BigRational::new_raw(1.into(), 0.into())),
187 /// Err(LaError::NonFinite {
188 /// location: NonFiniteLocation::MatrixCell { row: 0, col: 1, .. },
189 /// origin: NonFiniteOrigin::Input,
190 /// ..
191 /// })
192 /// );
193 /// assert_eq!(matrix, before);
194 /// # Ok(())
195 /// # }
196 /// ```
197 ///
198 /// # Errors
199 /// Returns [`LaError::IndexOutOfBounds`] when `(row, col)` lies outside the
200 /// matrix, or [`LaError::NonFinite`] when `value` has a raw zero
201 /// denominator.
202 pub fn set(&mut self, row: usize, col: usize, value: BigRational) -> Result<(), LaError> {
203 if row >= D || col >= D {
204 return Err(LaError::index_out_of_bounds(row, col, D));
205 }
206 if value.denom().sign() == Sign::NoSign {
207 return Err(LaError::non_finite_input_matrix(row, col));
208 }
209 self.rows[row][col] = canonicalize_rational(value);
210 Ok(())
211 }
212
213 /// Return the provably exact determinant sign.
214 ///
215 /// This path clears denominators and reads the sign of the resulting
216 /// integer determinant. It does not construct a rational determinant.
217 /// For D=0, the empty-product determinant has positive sign.
218 /// Use [`det`](Self::det) when the determinant value is also needed.
219 ///
220 /// # Examples
221 /// ```
222 /// use la_stack::prelude::*;
223 ///
224 /// # fn main() -> Result<(), LaError> {
225 /// let mut matrix = RationalMatrix::<2>::zero();
226 /// assert_eq!(matrix.det_sign(), DeterminantSign::Zero);
227 /// matrix.set(0, 1, BigRational::new(1.into(), 3.into()))?;
228 /// matrix.set(1, 0, BigRational::new(1.into(), 2.into()))?;
229 /// // The exact determinant is -1/6, so its sign is negative.
230 /// assert_eq!(matrix.det_sign(), DeterminantSign::Negative);
231 /// # Ok(())
232 /// # }
233 /// ```
234 pub fn det_sign(&self) -> DeterminantSign {
235 let (integer_rows, _) = self.integer_rows();
236 match det_big_int(integer_rows).sign() {
237 Sign::Minus => DeterminantSign::Negative,
238 Sign::NoSign => DeterminantSign::Zero,
239 Sign::Plus => DeterminantSign::Positive,
240 }
241 }
242
243 /// Return the exact determinant.
244 ///
245 /// Denominators are cleared independently per row. If row `i` uses
246 /// positive scale `sᵢ`, the integer determinant is divided by `∏ᵢ sᵢ`.
247 /// For D=0, this returns the empty-product determinant `1`.
248 /// Use [`det_sign`](Self::det_sign) when only the sign is needed, or
249 /// [`ExactF64Conversion`](crate::ExactF64Conversion) to convert this result
250 /// under an explicit strict or rounded binary64 contract.
251 ///
252 /// # Examples
253 /// ```
254 /// use la_stack::prelude::*;
255 ///
256 /// # fn main() -> Result<(), LaError> {
257 /// let mut matrix = RationalMatrix::<2>::zero();
258 /// matrix.set(0, 0, BigRational::new(1.into(), 3.into()))?;
259 /// matrix.set(1, 1, BigRational::from_integer(2.into()))?;
260 /// let determinant = matrix.det();
261 /// assert_eq!(determinant, BigRational::new(2.into(), 3.into()));
262 /// // Conversion rounds only after the exact determinant has been computed.
263 /// assert_eq!(determinant.to_rounded_f64()?, 2.0 / 3.0);
264 /// # Ok(())
265 /// # }
266 /// ```
267 #[must_use]
268 pub fn det(&self) -> BigRational {
269 let (integer_rows, row_scales) = self.integer_rows();
270 let determinant_denominator = row_scales.iter().product();
271 BigRational::new(det_big_int(integer_rows), determinant_denominator)
272 }
273
274 /// Solve `A x = b` exactly.
275 ///
276 /// Upper/lower triangular systems, including row permutations, use direct
277 /// exact substitution after `O(D²)` structure recognition. General systems
278 /// scale each augmented row `[A[row, :], rhs[row]]` to integers once, then
279 /// use fraction-free Bareiss forward elimination in [`BigInt`] followed by
280 /// `O(D²)` rational back-substitution. No binary64 conversion is performed.
281 /// For D=0, the empty matrix and vector have the unique empty solution.
282 ///
283 /// # Examples
284 /// ```
285 /// use la_stack::prelude::*;
286 ///
287 /// # fn main() -> Result<(), LaError> {
288 /// let zero = BigRational::from_integer(0.into());
289 /// let one = BigRational::from_integer(1.into());
290 /// let matrix = RationalMatrix::<2>::try_from_rows([
291 /// [BigRational::new(1.into(), 2.into()), zero.clone()],
292 /// [zero, BigRational::new(1.into(), 3.into())],
293 /// ])?;
294 /// let rhs = RationalVector::try_new([one.clone(), one])?;
295 ///
296 /// let solution = matrix.solve(&rhs)?.try_to_f64()?.into_array();
297 /// assert_eq!(solution, [2.0, 3.0]);
298 /// # Ok(())
299 /// # }
300 /// ```
301 ///
302 /// # Errors
303 /// Returns [`LaError::Singular`] with exact-singularity metadata when a
304 /// pivot column contains no non-zero entry.
305 pub fn solve(&self, rhs: &RationalVector<D>) -> Result<RationalVector<D>, LaError> {
306 if let Some(solution) = self.solve_triangular(rhs) {
307 return Ok(solution);
308 }
309 let row_scales: [BigInt; D] = from_fn(|row| {
310 common_denominator(
311 self.rows[row]
312 .iter()
313 .chain(core::iter::once(&rhs.data[row])),
314 )
315 });
316 let integer_rows =
317 from_fn(|row| from_fn(|col| integer_at_scale(&self.rows[row][col], &row_scales[row])));
318 let integer_rhs = from_fn(|row| integer_at_scale(&rhs.data[row], &row_scales[row]));
319 solve_big_int(integer_rows, integer_rhs).map(RationalVector::from_canonical_array)
320 }
321
322 /// A distinct first (upper) or last (lower) non-zero column in every row
323 /// proves a row permutation to a triangular matrix with non-zero diagonal.
324 ///
325 /// D distinct indices in 0..D cover every column; no matrix copy is needed.
326 /// Dense inputs reject at the second repeated edge column.
327 fn triangular_row_order(&self, upper: bool) -> Option<[usize; D]> {
328 let mut order = [0; D];
329 let mut occupied = [false; D];
330 for (row, entries) in self.rows.iter().enumerate() {
331 let column = if upper {
332 entries.iter().position(|entry| !entry.is_zero())
333 } else {
334 entries.iter().rposition(|entry| !entry.is_zero())
335 }?;
336 if occupied[column] {
337 return None;
338 }
339 occupied[column] = true;
340 order[column] = row;
341 }
342 Some(order)
343 }
344
345 /// Substitute only after every pivot has been proven non-zero.
346 ///
347 /// Accumulate each residual as an integer numerator and positive denominator,
348 /// using cross products instead of reducing every multiply and subtraction.
349 /// Reduce once after dividing by the proven non-zero pivot; stored solution
350 /// components are always canonical. Failed recognition returns to Bareiss,
351 /// preserving its singular pivot metadata.
352 fn solve_triangular(&self, rhs: &RationalVector<D>) -> Option<RationalVector<D>> {
353 let (order, upper) = self
354 .triangular_row_order(true)
355 .map(|order| (order, true))
356 .or_else(|| self.triangular_row_order(false).map(|order| (order, false)))?;
357 let mut solution: [BigRational; D] = from_fn(|_| BigRational::zero());
358 for step in 0..D {
359 let column = if upper { D - 1 - step } else { step };
360 let row = order[column];
361 let mut numerator = rhs.data[row].numer().clone();
362 let mut denominator = rhs.data[row].denom().clone();
363 let known = if upper { column + 1..D } else { 0..column };
364 for (coefficient, component) in
365 self.rows[row][known.clone()].iter().zip(&solution[known])
366 {
367 if coefficient.is_zero() || component.is_zero() {
368 continue;
369 }
370 let product_numerator = coefficient.numer() * component.numer();
371 let product_denominator = coefficient.denom() * component.denom();
372 if denominator == product_denominator {
373 numerator -= product_numerator;
374 } else {
375 numerator = numerator * &product_denominator - product_numerator * &denominator;
376 denominator *= product_denominator;
377 }
378 }
379 let pivot = &self.rows[row][column];
380 solution[column] =
381 BigRational::new(numerator * pivot.denom(), denominator * pivot.numer());
382 }
383 Some(RationalVector::from_canonical_array(solution))
384 }
385
386 /// Clear matrix denominators with one positive common denominator per row.
387 fn integer_rows(&self) -> ([[BigInt; D]; D], [BigInt; D]) {
388 let row_scales: [BigInt; D] = from_fn(|row| common_denominator(self.rows[row].iter()));
389 let integer_rows =
390 from_fn(|row| from_fn(|col| integer_at_scale(&self.rows[row][col], &row_scales[row])));
391 (integer_rows, row_scales)
392 }
393}
394
395impl<const D: usize> RationalVector<D> {
396 /// Wrap values produced by `BigRational` arithmetic, which preserves the
397 /// canonical representation established at public input boundaries.
398 pub(crate) const fn from_canonical_array(data: [BigRational; D]) -> Self {
399 Self { data }
400 }
401
402 /// Try to create an exact vector from rational storage.
403 ///
404 /// Raw, non-reduced values and negative denominators are accepted and
405 /// stored canonically. A raw zero denominator is rejected.
406 ///
407 /// # Examples
408 /// ```
409 /// use la_stack::prelude::*;
410 ///
411 /// # fn main() -> Result<(), LaError> {
412 /// let rhs = RationalVector::<2>::try_new([
413 /// BigRational::new(1.into(), 2.into()),
414 /// BigRational::from_integer(3.into()),
415 /// ])?;
416 /// assert_eq!(rhs.try_to_f64()?.into_array(), [0.5, 3.0]);
417 /// # Ok(())
418 /// # }
419 /// ```
420 ///
421 /// # Errors
422 /// Returns [`LaError::NonFinite`] at the first vector entry whose raw
423 /// rational denominator is zero.
424 pub fn try_new(data: [BigRational; D]) -> Result<Self, LaError> {
425 for (index, value) in data.iter().enumerate() {
426 if value.denom().sign() == Sign::NoSign {
427 return Err(LaError::non_finite_input_vector(index));
428 }
429 }
430 Ok(Self {
431 data: data.map(canonicalize_rational),
432 })
433 }
434
435 /// Try to create an exact vector by evaluating a function at every index.
436 ///
437 /// The function is evaluated once per entry in increasing index order.
438 /// All entries are generated before validation and canonicalization, as
439 /// in [`try_new`](Self::try_new).
440 ///
441 /// # Examples
442 /// ```
443 /// use la_stack::prelude::*;
444 ///
445 /// # fn main() -> Result<(), LaError> {
446 /// let numerators = [1, 2, 3];
447 /// let rhs = RationalVector::<3>::try_from_fn(|index| {
448 /// BigRational::new(numerators[index].into(), 2.into())
449 /// })?;
450 /// assert_eq!(rhs.try_to_f64()?.into_array(), [0.5, 1.0, 1.5]);
451 /// # Ok(())
452 /// # }
453 /// ```
454 ///
455 /// # Errors
456 /// Returns [`LaError::NonFinite`] at the first generated entry whose raw
457 /// rational denominator is zero.
458 pub fn try_from_fn(make_entry: impl FnMut(usize) -> BigRational) -> Result<Self, LaError> {
459 Self::try_new(from_fn(make_entry))
460 }
461
462 /// Return the all-zero exact vector.
463 pub fn zero() -> Self {
464 Self {
465 data: from_fn(|_| BigRational::from_integer(BigInt::from(0))),
466 }
467 }
468
469 /// Borrow the exact backing array.
470 #[must_use]
471 pub const fn as_array(&self) -> &[BigRational; D] {
472 &self.data
473 }
474
475 /// Consume the vector and return its exact backing array.
476 #[must_use]
477 pub fn into_array(self) -> [BigRational; D] {
478 self.data
479 }
480
481 /// Borrow one entry, returning `None` for an out-of-bounds index.
482 #[must_use]
483 pub fn get(&self, index: usize) -> Option<&BigRational> {
484 self.data.get(index)
485 }
486}
487
488/// Reduce one validated rational and make its denominator positive before
489/// publishing it through the exact-input storage types.
490fn canonicalize_rational(value: BigRational) -> BigRational {
491 let (numerator, denominator) = value.into_raw();
492 BigRational::new(numerator, denominator)
493}
494
495/// Return the least common multiple of canonical positive denominators.
496fn common_denominator<'a>(values: impl Iterator<Item = &'a BigRational>) -> BigInt {
497 values.fold(BigInt::from(1), |scale, value| {
498 least_common_multiple(scale, value.denom())
499 })
500}
501
502/// Convert a canonical rational to an integer using a positive divisible scale.
503///
504/// Matrix and vector construction already prove that the denominator is positive.
505fn integer_at_scale(value: &BigRational, scale: &BigInt) -> BigInt {
506 value.numer() * (scale / value.denom())
507}
508
509/// Return the positive least common multiple of two positive integers.
510///
511/// `lcm(a, b) = (a / gcd(a, b)) × b`; dividing before multiplying avoids
512/// forming the larger intermediate `a × b`. The resulting positive scale
513/// clears both denominators without changing determinant sign.
514fn least_common_multiple(lhs: BigInt, rhs: &BigInt) -> BigInt {
515 // Integer entries and repeated denominators leave the current scale unchanged.
516 if rhs.is_one() || lhs == *rhs {
517 return lhs;
518 }
519 // The first non-integer entry establishes the scale without a GCD.
520 if lhs.is_one() {
521 return rhs.clone();
522 }
523 let gcd = greatest_common_divisor(lhs.clone(), rhs.clone());
524 (lhs / gcd) * rhs
525}
526
527/// Euclidean greatest common divisor for positive integers.
528fn greatest_common_divisor(mut lhs: BigInt, mut rhs: BigInt) -> BigInt {
529 let zero = BigInt::from(0);
530 while rhs != zero {
531 let remainder = &lhs % &rhs;
532 lhs = rhs;
533 rhs = remainder;
534 }
535 lhs
536}
537
538#[cfg(test)]
539mod tests {
540 use pastey::paste;
541
542 use super::*;
543 use crate::{
544 ExactF64Conversion, NonFiniteLocation, NonFiniteOrigin, SingularityReason,
545 UnrepresentableReason,
546 };
547
548 fn ratio(numerator: i64, denominator: i64) -> BigRational {
549 BigRational::new(BigInt::from(numerator), BigInt::from(denominator))
550 }
551
552 #[test]
553 fn exact_input_preserves_non_binary64_coefficients() {
554 let matrix = RationalMatrix::<2>::try_from_rows([
555 [ratio(1, 3), ratio(1, 10)],
556 [ratio(2, 7), ratio(3, 11)],
557 ])
558 .unwrap();
559
560 assert_eq!(matrix.det(), ratio(24, 385));
561 assert_eq!(matrix.det_sign(), DeterminantSign::Positive);
562 }
563
564 #[test]
565 fn determinant_handles_pivoting_and_singularity() {
566 let pivoting = RationalMatrix::<3>::try_from_rows([
567 [ratio(0, 1), ratio(1, 2), ratio(0, 1)],
568 [ratio(2, 3), ratio(0, 1), ratio(0, 1)],
569 [ratio(0, 1), ratio(0, 1), ratio(3, 5)],
570 ])
571 .unwrap();
572 assert_eq!(pivoting.det(), ratio(-1, 5));
573 assert_eq!(pivoting.det_sign(), DeterminantSign::Negative);
574
575 let singular = RationalMatrix::<3>::try_from_rows([
576 [ratio(1, 2), ratio(1, 3), ratio(1, 5)],
577 [ratio(1, 1), ratio(2, 3), ratio(2, 5)],
578 [ratio(0, 1), ratio(1, 1), ratio(1, 1)],
579 ])
580 .unwrap();
581 assert_eq!(singular.det(), ratio(0, 1));
582 assert_eq!(singular.det_sign(), DeterminantSign::Zero);
583 }
584
585 #[test]
586 fn row_clearing_accepts_raw_signs_factors_and_dyadic_exponents() {
587 let matrix = RationalMatrix::<2>::try_from_rows([
588 [
589 BigRational::new_raw(BigInt::from(6), BigInt::from(-8)),
590 BigRational::new_raw(BigInt::from(3), BigInt::from(1_u8) << 80_u32),
591 ],
592 [
593 BigRational::new_raw(BigInt::from(-10), BigInt::from(-20)),
594 BigRational::new_raw(BigInt::from(14), BigInt::from(21)),
595 ],
596 ])
597 .unwrap();
598
599 let expected =
600 ratio(-1, 2) - BigRational::new(BigInt::from(3), BigInt::from(1_u8) << 81_u32);
601 assert_eq!(matrix.det(), expected);
602 assert_eq!(matrix.det_sign(), DeterminantSign::Negative);
603 assert_eq!(matrix.as_rows()[0][0], ratio(-3, 4));
604 assert_eq!(matrix.as_rows()[1][0], ratio(1, 2));
605 assert_eq!(matrix.as_rows()[1][1], ratio(2, 3));
606 }
607
608 #[test]
609 fn construction_and_set_store_equal_quotients_identically() {
610 let huge_factor = BigInt::from(1_u8) << 256_u32;
611 let raw = RationalMatrix::<2>::try_from_rows([
612 [
613 BigRational::new_raw(huge_factor.clone(), &huge_factor * 2_u8),
614 BigRational::new_raw(BigInt::from(0), -&huge_factor),
615 ],
616 [ratio(0, 1), ratio(1, 1)],
617 ])
618 .unwrap();
619 let canonical = RationalMatrix::<2>::try_from_rows([
620 [ratio(1, 2), ratio(0, 1)],
621 [ratio(0, 1), ratio(1, 1)],
622 ])
623 .unwrap();
624
625 assert_eq!(raw, canonical);
626 assert_eq!(raw.as_rows()[0][1].denom(), &BigInt::from(1));
627
628 let mut updated = RationalMatrix::<2>::zero();
629 updated
630 .set(
631 0,
632 0,
633 BigRational::new_raw(&huge_factor * 3_u8, -&huge_factor * 6_u8),
634 )
635 .unwrap();
636 assert_eq!(updated.get(0, 0), Some(&ratio(-1, 2)));
637 }
638
639 #[test]
640 fn exact_solve_scales_matrix_and_rhs_together() {
641 let matrix = RationalMatrix::<3>::try_from_rows([
642 [ratio(0, 1), ratio(1, 3), ratio(1, 5)],
643 [ratio(2, 7), ratio(1, 11), ratio(0, 1)],
644 [ratio(1, 13), ratio(0, 1), ratio(3, 2)],
645 ])
646 .unwrap();
647 let expected = [ratio(2, 3), ratio(-4, 5), ratio(7, 9)];
648 let rhs = RationalVector::try_from_fn(|row| {
649 matrix.as_rows()[row]
650 .iter()
651 .zip(expected.iter())
652 .map(|(coefficient, solution)| coefficient * solution)
653 .sum()
654 })
655 .unwrap();
656
657 assert_eq!(matrix.solve(&rhs).unwrap().into_array(), expected);
658 }
659
660 #[test]
661 fn exact_solve_accepts_signed_unreduced_matrix_and_rhs_values() {
662 let matrix = RationalMatrix::<2>::try_from_rows([
663 [
664 BigRational::new_raw(BigInt::from(10), BigInt::from(20)),
665 BigRational::new_raw(BigInt::from(-5), BigInt::from(-15)),
666 ],
667 [
668 BigRational::new_raw(BigInt::from(-6), BigInt::from(-15)),
669 BigRational::new_raw(BigInt::from(9), BigInt::from(21)),
670 ],
671 ])
672 .unwrap();
673 let rhs = RationalVector::try_new([
674 BigRational::new_raw(BigInt::from(0), BigInt::from(-99)),
675 BigRational::new_raw(BigInt::from(34), BigInt::from(-70)),
676 ])
677 .unwrap();
678
679 let solution = matrix.solve(&rhs).unwrap();
680 assert_eq!(solution.as_array(), &[ratio(2, 1), ratio(-3, 1)]);
681 assert_eq!(
682 &matrix.as_rows()[0][0] * &solution.as_array()[0]
683 + &matrix.as_rows()[0][1] * &solution.as_array()[1],
684 rhs.as_array()[0]
685 );
686 assert_eq!(
687 &matrix.as_rows()[1][0] * &solution.as_array()[0]
688 + &matrix.as_rows()[1][1] * &solution.as_array()[1],
689 rhs.as_array()[1]
690 );
691 }
692
693 #[test]
694 fn singular_solve_preserves_exact_pivot_metadata() {
695 let matrix = RationalMatrix::<2>::try_from_rows([
696 [ratio(1, 2), ratio(1, 3)],
697 [ratio(1, 1), ratio(2, 3)],
698 ])
699 .unwrap();
700 let rhs = RationalVector::<2>::zero();
701
702 assert!(matches!(
703 matrix.solve(&rhs),
704 Err(LaError::Singular {
705 pivot_col: 1,
706 reason: SingularityReason::Exact,
707 ..
708 })
709 ));
710 }
711
712 #[test]
713 fn constructors_reject_raw_zero_denominators_with_locations() {
714 let matrix_error = RationalMatrix::<1>::try_from_rows([[BigRational::new_raw(
715 BigInt::from(1),
716 BigInt::from(0),
717 )]])
718 .unwrap_err();
719 assert!(matches!(
720 matrix_error,
721 LaError::NonFinite {
722 location: NonFiniteLocation::MatrixCell { row: 0, col: 0, .. },
723 origin: NonFiniteOrigin::Input,
724 ..
725 }
726 ));
727
728 let vector_error =
729 RationalVector::<1>::try_new([BigRational::new_raw(BigInt::from(1), BigInt::from(0))])
730 .unwrap_err();
731 assert!(matches!(
732 vector_error,
733 LaError::NonFinite {
734 location: NonFiniteLocation::VectorEntry { index: 0, .. },
735 origin: NonFiniteOrigin::Input,
736 ..
737 }
738 ));
739 }
740
741 macro_rules! gen_rejected_set_is_failure_atomic_tests {
742 ($d:literal) => {
743 paste! {
744 #[test]
745 fn [<rejected_set_is_failure_atomic_ $d d>]() {
746 let mut matrix = RationalMatrix::<$d>::try_from_fn(|row, col| {
747 BigRational::from_integer(BigInt::from(row * $d + col + 1))
748 })
749 .unwrap();
750 let original = matrix.clone();
751
752 assert!(matches!(
753 matrix.set($d, 0, ratio(1, 2)),
754 Err(LaError::IndexOutOfBounds {
755 row: $d,
756 col: 0,
757 dim: $d,
758 ..
759 })
760 ));
761 assert_eq!(matrix, original);
762
763 let zero_denominator =
764 BigRational::new_raw(BigInt::from(1), BigInt::from(0));
765 assert!(matches!(
766 matrix.set(0, 1, zero_denominator),
767 Err(LaError::NonFinite {
768 location: NonFiniteLocation::MatrixCell { row: 0, col: 1, .. },
769 origin: NonFiniteOrigin::Input,
770 ..
771 })
772 ));
773 assert_eq!(matrix, original);
774 }
775 }
776 };
777 }
778
779 gen_rejected_set_is_failure_atomic_tests!(2);
780 gen_rejected_set_is_failure_atomic_tests!(3);
781 gen_rejected_set_is_failure_atomic_tests!(4);
782 gen_rejected_set_is_failure_atomic_tests!(5);
783
784 #[test]
785 fn zero_dimension_uses_empty_determinant_and_unique_solve() {
786 let matrix = RationalMatrix::<0>::zero();
787 assert_eq!(matrix.det(), ratio(1, 1));
788 assert_eq!(matrix.det_sign(), DeterminantSign::Positive);
789 assert_eq!(
790 matrix.solve(&RationalVector::zero()).unwrap().into_array(),
791 []
792 );
793 }
794
795 #[test]
796 fn rational_vector_conversion_is_explicit() {
797 let vector = RationalVector::<2>::try_new([ratio(1, 2), ratio(1, 3)]).unwrap();
798 assert!(matches!(
799 vector.try_to_f64(),
800 Err(LaError::Unrepresentable {
801 index: Some(1),
802 reason: UnrepresentableReason::RequiresRounding,
803 ..
804 })
805 ));
806 let rounded = vector.to_rounded_f64().unwrap().into_array();
807 assert_eq!(rounded[0].to_bits(), 0.5_f64.to_bits());
808 assert_eq!(rounded[1].to_bits(), (1.0_f64 / 3.0).to_bits());
809 }
810
811 macro_rules! gen_pivoting_tests {
812 ($d:literal) => {
813 paste! {
814 #[test]
815 fn [<pivoting_determinant_and_solve_ $d d>]() {
816 let matrix = RationalMatrix::<$d>::try_from_fn(|row, col| {
817 let is_permutation_entry = (row == 0 && col == 1)
818 || (row == 1 && col == 0)
819 || (row >= 2 && row == col);
820 BigRational::from_integer(BigInt::from(u8::from(is_permutation_entry)))
821 })
822 .unwrap();
823 let expected = from_fn(|index| {
824 BigRational::new(BigInt::from(1), BigInt::from(index + 2))
825 });
826 let rhs = RationalVector::try_from_fn(|row| {
827 matrix.as_rows()[row]
828 .iter()
829 .zip(expected.iter())
830 .map(|(coefficient, component)| coefficient * component)
831 .sum()
832 })
833 .unwrap();
834
835 assert_eq!(matrix.det_sign(), DeterminantSign::Negative);
836 assert_eq!(matrix.det(), ratio(-1, 1));
837 assert_eq!(matrix.solve(&rhs).unwrap().into_array(), expected);
838 }
839 }
840 };
841 }
842
843 gen_pivoting_tests!(2);
844 gen_pivoting_tests!(3);
845 gen_pivoting_tests!(4);
846 gen_pivoting_tests!(5);
847 gen_pivoting_tests!(6);
848 gen_pivoting_tests!(7);
849 gen_pivoting_tests!(8);
850}