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