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²)` integer back-substitution over a shared denominator, reducing once
281 /// per solved row and canonicalizing final outputs. No binary64 conversion
282 /// is performed.
283 /// For D=0, the empty matrix and vector have the unique empty solution.
284 ///
285 /// # Examples
286 /// ```
287 /// use la_stack::prelude::*;
288 ///
289 /// # fn main() -> Result<(), LaError> {
290 /// let zero = BigRational::from_integer(0.into());
291 /// let one = BigRational::from_integer(1.into());
292 /// let matrix = RationalMatrix::<2>::try_from_rows([
293 /// [BigRational::new(1.into(), 2.into()), zero.clone()],
294 /// [zero, BigRational::new(1.into(), 3.into())],
295 /// ])?;
296 /// let rhs = RationalVector::try_new([one.clone(), one])?;
297 ///
298 /// let solution = matrix.solve(&rhs)?.try_to_f64()?.into_array();
299 /// assert_eq!(solution, [2.0, 3.0]);
300 /// # Ok(())
301 /// # }
302 /// ```
303 ///
304 /// # Errors
305 /// Returns [`LaError::Singular`] with exact-singularity metadata when a
306 /// pivot column contains no non-zero entry.
307 pub fn solve(&self, rhs: &RationalVector<D>) -> Result<RationalVector<D>, LaError> {
308 if let Some(solution) = self.solve_triangular(rhs) {
309 return Ok(solution);
310 }
311 let row_scales: [BigInt; D] = from_fn(|row| {
312 common_denominator(
313 self.rows[row]
314 .iter()
315 .chain(core::iter::once(&rhs.data[row])),
316 )
317 });
318 let integer_rows =
319 from_fn(|row| from_fn(|col| integer_at_scale(&self.rows[row][col], &row_scales[row])));
320 let integer_rhs = from_fn(|row| integer_at_scale(&rhs.data[row], &row_scales[row]));
321 solve_big_int(integer_rows, integer_rhs).map(RationalVector::from_canonical_array)
322 }
323
324 /// A distinct first (upper) or last (lower) non-zero column in every row
325 /// proves a row permutation to a triangular matrix with non-zero diagonal.
326 ///
327 /// D distinct indices in 0..D cover every column; no matrix copy is needed.
328 /// Dense inputs reject at the second repeated edge column.
329 fn triangular_row_order(&self, upper: bool) -> Option<[usize; D]> {
330 let mut order = [0; D];
331 let mut occupied = [false; D];
332 for (row, entries) in self.rows.iter().enumerate() {
333 let column = if upper {
334 entries.iter().position(|entry| !entry.is_zero())
335 } else {
336 entries.iter().rposition(|entry| !entry.is_zero())
337 }?;
338 if occupied[column] {
339 return None;
340 }
341 occupied[column] = true;
342 order[column] = row;
343 }
344 Some(order)
345 }
346
347 /// Substitute only after every pivot has been proven non-zero.
348 ///
349 /// Accumulate each residual as an integer numerator and positive denominator,
350 /// using cross products instead of reducing every multiply and subtraction.
351 /// Reduce once after dividing by the proven non-zero pivot; stored solution
352 /// components are always canonical. Failed recognition returns to Bareiss,
353 /// preserving its singular pivot metadata.
354 fn solve_triangular(&self, rhs: &RationalVector<D>) -> Option<RationalVector<D>> {
355 let (order, upper) = self
356 .triangular_row_order(true)
357 .map(|order| (order, true))
358 .or_else(|| self.triangular_row_order(false).map(|order| (order, false)))?;
359 let mut solution: [BigRational; D] = from_fn(|_| BigRational::zero());
360 for step in 0..D {
361 let column = if upper { D - 1 - step } else { step };
362 let row = order[column];
363 let mut numerator = rhs.data[row].numer().clone();
364 let mut denominator = rhs.data[row].denom().clone();
365 let known = if upper { column + 1..D } else { 0..column };
366 for (coefficient, component) in
367 self.rows[row][known.clone()].iter().zip(&solution[known])
368 {
369 if coefficient.is_zero() || component.is_zero() {
370 continue;
371 }
372 let product_numerator = coefficient.numer() * component.numer();
373 let product_denominator = coefficient.denom() * component.denom();
374 if denominator == product_denominator {
375 numerator -= product_numerator;
376 } else {
377 numerator = numerator * &product_denominator - product_numerator * &denominator;
378 denominator *= product_denominator;
379 }
380 }
381 let pivot = &self.rows[row][column];
382 solution[column] =
383 BigRational::new(numerator * pivot.denom(), denominator * pivot.numer());
384 }
385 Some(RationalVector::from_canonical_array(solution))
386 }
387
388 /// Clear matrix denominators with one positive common denominator per row.
389 fn integer_rows(&self) -> ([[BigInt; D]; D], [BigInt; D]) {
390 let row_scales: [BigInt; D] = from_fn(|row| common_denominator(self.rows[row].iter()));
391 let integer_rows =
392 from_fn(|row| from_fn(|col| integer_at_scale(&self.rows[row][col], &row_scales[row])));
393 (integer_rows, row_scales)
394 }
395}
396
397impl<const D: usize> RationalVector<D> {
398 /// Wrap values produced by `BigRational` arithmetic, which preserves the
399 /// canonical representation established at public input boundaries.
400 pub(crate) const fn from_canonical_array(data: [BigRational; D]) -> Self {
401 Self { data }
402 }
403
404 /// Try to create an exact vector from rational storage.
405 ///
406 /// Raw, non-reduced values and negative denominators are accepted and
407 /// stored canonically. A raw zero denominator is rejected.
408 ///
409 /// # Examples
410 /// ```
411 /// use la_stack::prelude::*;
412 ///
413 /// # fn main() -> Result<(), LaError> {
414 /// let rhs = RationalVector::<2>::try_new([
415 /// BigRational::new(1.into(), 2.into()),
416 /// BigRational::from_integer(3.into()),
417 /// ])?;
418 /// assert_eq!(rhs.try_to_f64()?.into_array(), [0.5, 3.0]);
419 /// # Ok(())
420 /// # }
421 /// ```
422 ///
423 /// # Errors
424 /// Returns [`LaError::NonFinite`] at the first vector entry whose raw
425 /// rational denominator is zero.
426 pub fn try_new(data: [BigRational; D]) -> Result<Self, LaError> {
427 for (index, value) in data.iter().enumerate() {
428 if value.denom().sign() == Sign::NoSign {
429 return Err(LaError::non_finite_input_vector(index));
430 }
431 }
432 Ok(Self {
433 data: data.map(canonicalize_rational),
434 })
435 }
436
437 /// Try to create an exact vector by evaluating a function at every index.
438 ///
439 /// The function is evaluated once per entry in increasing index order.
440 /// All entries are generated before validation and canonicalization, as
441 /// in [`try_new`](Self::try_new).
442 ///
443 /// # Examples
444 /// ```
445 /// use la_stack::prelude::*;
446 ///
447 /// # fn main() -> Result<(), LaError> {
448 /// let numerators = [1, 2, 3];
449 /// let rhs = RationalVector::<3>::try_from_fn(|index| {
450 /// BigRational::new(numerators[index].into(), 2.into())
451 /// })?;
452 /// assert_eq!(rhs.try_to_f64()?.into_array(), [0.5, 1.0, 1.5]);
453 /// # Ok(())
454 /// # }
455 /// ```
456 ///
457 /// # Errors
458 /// Returns [`LaError::NonFinite`] at the first generated entry whose raw
459 /// rational denominator is zero.
460 pub fn try_from_fn(make_entry: impl FnMut(usize) -> BigRational) -> Result<Self, LaError> {
461 Self::try_new(from_fn(make_entry))
462 }
463
464 /// Return the all-zero exact vector.
465 pub fn zero() -> Self {
466 Self {
467 data: from_fn(|_| BigRational::from_integer(BigInt::from(0))),
468 }
469 }
470
471 /// Borrow the exact backing array.
472 #[must_use]
473 pub const fn as_array(&self) -> &[BigRational; D] {
474 &self.data
475 }
476
477 /// Consume the vector and return its exact backing array.
478 #[must_use]
479 pub fn into_array(self) -> [BigRational; D] {
480 self.data
481 }
482
483 /// Borrow one entry, returning `None` for an out-of-bounds index.
484 #[must_use]
485 pub fn get(&self, index: usize) -> Option<&BigRational> {
486 self.data.get(index)
487 }
488}
489
490/// Reduce one validated rational and make its denominator positive before
491/// publishing it through the exact-input storage types.
492fn canonicalize_rational(value: BigRational) -> BigRational {
493 let (numerator, denominator) = value.into_raw();
494 BigRational::new(numerator, denominator)
495}
496
497/// Return the least common multiple of canonical positive denominators.
498fn common_denominator<'a>(values: impl Iterator<Item = &'a BigRational>) -> BigInt {
499 values.fold(BigInt::from(1), |scale, value| {
500 least_common_multiple(scale, value.denom())
501 })
502}
503
504/// Convert a canonical rational to an integer using a positive divisible scale.
505///
506/// Matrix and vector construction already prove that the denominator is positive.
507fn integer_at_scale(value: &BigRational, scale: &BigInt) -> BigInt {
508 value.numer() * (scale / value.denom())
509}
510
511/// Return the positive least common multiple of two positive integers.
512///
513/// `lcm(a, b) = (a / gcd(a, b)) × b`; dividing before multiplying avoids
514/// forming the larger intermediate `a × b`. The resulting positive scale
515/// clears both denominators without changing determinant sign.
516fn least_common_multiple(lhs: BigInt, rhs: &BigInt) -> BigInt {
517 // Integer entries and repeated denominators leave the current scale unchanged.
518 if rhs.is_one() || lhs == *rhs {
519 return lhs;
520 }
521 // The first non-integer entry establishes the scale without a GCD.
522 if lhs.is_one() {
523 return rhs.clone();
524 }
525 let gcd = greatest_common_divisor(lhs.clone(), rhs.clone());
526 (lhs / gcd) * rhs
527}
528
529/// Euclidean greatest common divisor for positive integers.
530fn greatest_common_divisor(mut lhs: BigInt, mut rhs: BigInt) -> BigInt {
531 let zero = BigInt::from(0);
532 while rhs != zero {
533 let remainder = &lhs % &rhs;
534 lhs = rhs;
535 rhs = remainder;
536 }
537 lhs
538}
539
540#[cfg(test)]
541mod tests {
542 use pastey::paste;
543
544 use super::*;
545 use crate::{
546 ExactF64Conversion, NonFiniteLocation, NonFiniteOrigin, SingularityReason,
547 UnrepresentableReason,
548 };
549
550 fn ratio(numerator: i64, denominator: i64) -> BigRational {
551 BigRational::new(BigInt::from(numerator), BigInt::from(denominator))
552 }
553
554 #[test]
555 fn exact_input_preserves_non_binary64_coefficients() {
556 let matrix = RationalMatrix::<2>::try_from_rows([
557 [ratio(1, 3), ratio(1, 10)],
558 [ratio(2, 7), ratio(3, 11)],
559 ])
560 .unwrap();
561
562 assert_eq!(matrix.det(), ratio(24, 385));
563 assert_eq!(matrix.det_sign(), DeterminantSign::Positive);
564 }
565
566 #[test]
567 fn determinant_handles_pivoting_and_singularity() {
568 let pivoting = RationalMatrix::<3>::try_from_rows([
569 [ratio(0, 1), ratio(1, 2), ratio(0, 1)],
570 [ratio(2, 3), ratio(0, 1), ratio(0, 1)],
571 [ratio(0, 1), ratio(0, 1), ratio(3, 5)],
572 ])
573 .unwrap();
574 assert_eq!(pivoting.det(), ratio(-1, 5));
575 assert_eq!(pivoting.det_sign(), DeterminantSign::Negative);
576
577 let singular = RationalMatrix::<3>::try_from_rows([
578 [ratio(1, 2), ratio(1, 3), ratio(1, 5)],
579 [ratio(1, 1), ratio(2, 3), ratio(2, 5)],
580 [ratio(0, 1), ratio(1, 1), ratio(1, 1)],
581 ])
582 .unwrap();
583 assert_eq!(singular.det(), ratio(0, 1));
584 assert_eq!(singular.det_sign(), DeterminantSign::Zero);
585 }
586
587 #[test]
588 fn row_clearing_accepts_raw_signs_factors_and_dyadic_exponents() {
589 let matrix = RationalMatrix::<2>::try_from_rows([
590 [
591 BigRational::new_raw(BigInt::from(6), BigInt::from(-8)),
592 BigRational::new_raw(BigInt::from(3), BigInt::from(1_u8) << 80_u32),
593 ],
594 [
595 BigRational::new_raw(BigInt::from(-10), BigInt::from(-20)),
596 BigRational::new_raw(BigInt::from(14), BigInt::from(21)),
597 ],
598 ])
599 .unwrap();
600
601 let expected =
602 ratio(-1, 2) - BigRational::new(BigInt::from(3), BigInt::from(1_u8) << 81_u32);
603 assert_eq!(matrix.det(), expected);
604 assert_eq!(matrix.det_sign(), DeterminantSign::Negative);
605 assert_eq!(matrix.as_rows()[0][0], ratio(-3, 4));
606 assert_eq!(matrix.as_rows()[1][0], ratio(1, 2));
607 assert_eq!(matrix.as_rows()[1][1], ratio(2, 3));
608 }
609
610 #[test]
611 fn construction_and_set_store_equal_quotients_identically() {
612 let huge_factor = BigInt::from(1_u8) << 256_u32;
613 let raw = RationalMatrix::<2>::try_from_rows([
614 [
615 BigRational::new_raw(huge_factor.clone(), &huge_factor * 2_u8),
616 BigRational::new_raw(BigInt::from(0), -&huge_factor),
617 ],
618 [ratio(0, 1), ratio(1, 1)],
619 ])
620 .unwrap();
621 let canonical = RationalMatrix::<2>::try_from_rows([
622 [ratio(1, 2), ratio(0, 1)],
623 [ratio(0, 1), ratio(1, 1)],
624 ])
625 .unwrap();
626
627 assert_eq!(raw, canonical);
628 assert_eq!(raw.as_rows()[0][1].denom(), &BigInt::from(1));
629
630 let mut updated = RationalMatrix::<2>::zero();
631 updated
632 .set(
633 0,
634 0,
635 BigRational::new_raw(&huge_factor * 3_u8, -&huge_factor * 6_u8),
636 )
637 .unwrap();
638 assert_eq!(updated.get(0, 0), Some(&ratio(-1, 2)));
639 }
640
641 #[test]
642 fn exact_solve_scales_matrix_and_rhs_together() {
643 let matrix = RationalMatrix::<3>::try_from_rows([
644 [ratio(0, 1), ratio(1, 3), ratio(1, 5)],
645 [ratio(2, 7), ratio(1, 11), ratio(0, 1)],
646 [ratio(1, 13), ratio(0, 1), ratio(3, 2)],
647 ])
648 .unwrap();
649 let expected = [ratio(2, 3), ratio(-4, 5), ratio(7, 9)];
650 let rhs = RationalVector::try_from_fn(|row| {
651 matrix.as_rows()[row]
652 .iter()
653 .zip(expected.iter())
654 .map(|(coefficient, solution)| coefficient * solution)
655 .sum()
656 })
657 .unwrap();
658
659 assert_eq!(matrix.solve(&rhs).unwrap().into_array(), expected);
660 }
661
662 #[test]
663 fn exact_solve_accepts_signed_unreduced_matrix_and_rhs_values() {
664 let matrix = RationalMatrix::<2>::try_from_rows([
665 [
666 BigRational::new_raw(BigInt::from(10), BigInt::from(20)),
667 BigRational::new_raw(BigInt::from(-5), BigInt::from(-15)),
668 ],
669 [
670 BigRational::new_raw(BigInt::from(-6), BigInt::from(-15)),
671 BigRational::new_raw(BigInt::from(9), BigInt::from(21)),
672 ],
673 ])
674 .unwrap();
675 let rhs = RationalVector::try_new([
676 BigRational::new_raw(BigInt::from(0), BigInt::from(-99)),
677 BigRational::new_raw(BigInt::from(34), BigInt::from(-70)),
678 ])
679 .unwrap();
680
681 let solution = matrix.solve(&rhs).unwrap();
682 assert_eq!(solution.as_array(), &[ratio(2, 1), ratio(-3, 1)]);
683 assert_eq!(
684 &matrix.as_rows()[0][0] * &solution.as_array()[0]
685 + &matrix.as_rows()[0][1] * &solution.as_array()[1],
686 rhs.as_array()[0]
687 );
688 assert_eq!(
689 &matrix.as_rows()[1][0] * &solution.as_array()[0]
690 + &matrix.as_rows()[1][1] * &solution.as_array()[1],
691 rhs.as_array()[1]
692 );
693 }
694
695 #[test]
696 fn singular_solve_preserves_exact_pivot_metadata() {
697 let matrix = RationalMatrix::<2>::try_from_rows([
698 [ratio(1, 2), ratio(1, 3)],
699 [ratio(1, 1), ratio(2, 3)],
700 ])
701 .unwrap();
702 let rhs = RationalVector::<2>::zero();
703
704 assert!(matches!(
705 matrix.solve(&rhs),
706 Err(LaError::Singular {
707 pivot_col: 1,
708 reason: SingularityReason::Exact,
709 ..
710 })
711 ));
712 }
713
714 #[test]
715 fn constructors_reject_raw_zero_denominators_with_locations() {
716 let matrix_error = RationalMatrix::<1>::try_from_rows([[BigRational::new_raw(
717 BigInt::from(1),
718 BigInt::from(0),
719 )]])
720 .unwrap_err();
721 assert!(matches!(
722 matrix_error,
723 LaError::NonFinite {
724 location: NonFiniteLocation::MatrixCell { row: 0, col: 0, .. },
725 origin: NonFiniteOrigin::Input,
726 ..
727 }
728 ));
729
730 let vector_error =
731 RationalVector::<1>::try_new([BigRational::new_raw(BigInt::from(1), BigInt::from(0))])
732 .unwrap_err();
733 assert!(matches!(
734 vector_error,
735 LaError::NonFinite {
736 location: NonFiniteLocation::VectorEntry { index: 0, .. },
737 origin: NonFiniteOrigin::Input,
738 ..
739 }
740 ));
741 }
742
743 macro_rules! gen_rejected_set_is_failure_atomic_tests {
744 ($d:literal) => {
745 paste! {
746 #[test]
747 fn [<rejected_set_is_failure_atomic_ $d d>]() {
748 let mut matrix = RationalMatrix::<$d>::try_from_fn(|row, col| {
749 BigRational::from_integer(BigInt::from(row * $d + col + 1))
750 })
751 .unwrap();
752 let original = matrix.clone();
753
754 assert!(matches!(
755 matrix.set($d, 0, ratio(1, 2)),
756 Err(LaError::IndexOutOfBounds {
757 row: $d,
758 col: 0,
759 dim: $d,
760 ..
761 })
762 ));
763 assert_eq!(matrix, original);
764
765 let zero_denominator =
766 BigRational::new_raw(BigInt::from(1), BigInt::from(0));
767 assert!(matches!(
768 matrix.set(0, 1, zero_denominator),
769 Err(LaError::NonFinite {
770 location: NonFiniteLocation::MatrixCell { row: 0, col: 1, .. },
771 origin: NonFiniteOrigin::Input,
772 ..
773 })
774 ));
775 assert_eq!(matrix, original);
776 }
777 }
778 };
779 }
780
781 gen_rejected_set_is_failure_atomic_tests!(2);
782 gen_rejected_set_is_failure_atomic_tests!(3);
783 gen_rejected_set_is_failure_atomic_tests!(4);
784 gen_rejected_set_is_failure_atomic_tests!(5);
785
786 #[test]
787 fn zero_dimension_uses_empty_determinant_and_unique_solve() {
788 let matrix = RationalMatrix::<0>::zero();
789 assert_eq!(matrix.det(), ratio(1, 1));
790 assert_eq!(matrix.det_sign(), DeterminantSign::Positive);
791 assert_eq!(
792 matrix.solve(&RationalVector::zero()).unwrap().into_array(),
793 []
794 );
795 }
796
797 #[test]
798 fn rational_vector_conversion_is_explicit() {
799 let vector = RationalVector::<2>::try_new([ratio(1, 2), ratio(1, 3)]).unwrap();
800 assert!(matches!(
801 vector.try_to_f64(),
802 Err(LaError::Unrepresentable {
803 index: Some(1),
804 reason: UnrepresentableReason::RequiresRounding,
805 ..
806 })
807 ));
808 let rounded = vector.to_rounded_f64().unwrap().into_array();
809 assert_eq!(rounded[0].to_bits(), 0.5_f64.to_bits());
810 assert_eq!(rounded[1].to_bits(), (1.0_f64 / 3.0).to_bits());
811 }
812
813 macro_rules! gen_pivoting_tests {
814 ($d:literal) => {
815 paste! {
816 #[test]
817 fn [<pivoting_determinant_and_solve_ $d d>]() {
818 let matrix = RationalMatrix::<$d>::try_from_fn(|row, col| {
819 let is_permutation_entry = (row == 0 && col == 1)
820 || (row == 1 && col == 0)
821 || (row >= 2 && row == col);
822 BigRational::from_integer(BigInt::from(u8::from(is_permutation_entry)))
823 })
824 .unwrap();
825 let expected = from_fn(|index| {
826 BigRational::new(BigInt::from(1), BigInt::from(index + 2))
827 });
828 let rhs = RationalVector::try_from_fn(|row| {
829 matrix.as_rows()[row]
830 .iter()
831 .zip(expected.iter())
832 .map(|(coefficient, component)| coefficient * component)
833 .sum()
834 })
835 .unwrap();
836
837 assert_eq!(matrix.det_sign(), DeterminantSign::Negative);
838 assert_eq!(matrix.det(), ratio(-1, 1));
839 assert_eq!(matrix.solve(&rhs).unwrap().into_array(), expected);
840 }
841 }
842 };
843 }
844
845 gen_pivoting_tests!(2);
846 gen_pivoting_tests!(3);
847 gen_pivoting_tests!(4);
848 gen_pivoting_tests!(5);
849 gen_pivoting_tests!(6);
850 gen_pivoting_tests!(7);
851 gen_pivoting_tests!(8);
852}