feanor_math/algorithms/linsolve/
mod.rs

1use std::alloc::{Allocator, Global};
2
3use crate::divisibility::DivisibilityRing;
4use crate::matrix::{AsPointerToSlice, SubmatrixMut};
5use crate::pid::PrincipalIdealRing;
6use crate::ring::*;
7use crate::rings::extension::extension_impl::FreeAlgebraImplBase;
8use crate::seq::VectorView;
9
10use super::convolution::ConvolutionAlgorithm;
11
12///
13/// Contains the algorithm to compute the "pre-smith" form of a matrix, which
14/// is sufficient for solving linear systems. Works over all principal ideal rings.
15/// 
16pub mod smith;
17///
18/// Contains the algorithm for solving linear systems over free ring extensions.
19/// 
20pub mod extension;
21///
22/// Contains algorithms related to Gaussian elimination. Note that standard functionality
23/// relating to solving linear systems is provided by [`LinSolveRing`] instead.
24/// 
25pub mod gauss;
26
27///
28/// Result of trying to solve a linear system.
29/// 
30/// Possible values are:
31///  - [`SolveResult::FoundUniqueSolution`]: The system is guaranteed to have a unique solution.
32///  - [`SolveResult::FoundSomeSolution`]: The system has at least one solution. This is also an allowed
33///    value for systems with a unique solution.
34///  - [`SolveResult::NoSolution`]: The system is unsolvable.
35/// 
36#[derive(Debug, Clone, PartialEq, Eq, Hash)]
37pub enum SolveResult {
38    /// The system has at least one solution. This is also an allowed
39    /// value for systems with a unique solution.
40    FoundSomeSolution,
41    /// The system is guaranteed to have a unique solution 
42    FoundUniqueSolution,
43    /// The system has no solution. In particular, the output matrix
44    /// is in an unspecified (but safe) state.
45    NoSolution
46}
47
48impl SolveResult {
49
50    ///
51    /// Returns whether some solution to the system has been found.
52    /// 
53    pub fn is_solved(&self) -> bool {
54        match self {
55            Self::FoundSomeSolution | Self::FoundUniqueSolution => true,
56            Self::NoSolution => false
57        }
58    }
59    
60    ///
61    /// Panics if the system does not have a solution.
62    /// 
63    pub fn assert_solved(&self) {
64        assert!(self.is_solved());
65    }
66}
67
68///
69/// Class for rings over which we can solve linear systems.
70/// 
71pub trait LinSolveRing: DivisibilityRing {
72
73    ///
74    /// Tries to find a matrix `X` such that `lhs * X = rhs`.
75    /// 
76    /// If a solution exists, it will be written to `out`. Otherwise, `out` will have an
77    /// unspecified (but valid) value after the function returns. Similarly, `lhs` and `rhs`
78    /// will be modified in an unspecified way, but its entries will always be valid ring elements.
79    /// 
80    /// Note that if a solution is found, either [`SolveResult::FoundSomeSolution`] or [`SolveResult::FoundUniqueSolution`]
81    /// are returned. If there are multiple solutions, only former will be returned. However, implementations
82    /// are free to also return [`SolveResult::FoundSomeSolution`] in cases where there is a unique solution.
83    /// 
84    fn solve_right<V1, V2, V3, A>(&self, lhs: SubmatrixMut<V1, Self::Element>, rhs: SubmatrixMut<V2, Self::Element>, out: SubmatrixMut<V3, Self::Element>, allocator: A) -> SolveResult
85        where V1: AsPointerToSlice<Self::Element>,
86            V2: AsPointerToSlice<Self::Element>,
87            V3: AsPointerToSlice<Self::Element>,
88            A: Allocator;
89}
90
91///
92/// [`RingStore`] corresponding to [`LinSolveRing`].
93/// 
94pub trait LinSolveRingStore: RingStore
95    where Self::Type: LinSolveRing
96{
97    ///
98    /// Solves a linear system `lhs * X = rhs`.
99    /// 
100    /// For details, see [`LinSolveRing::solve_right()`].
101    /// 
102    fn solve_right<V1, V2, V3>(&self, lhs: SubmatrixMut<V1, El<Self>>, rhs: SubmatrixMut<V2, El<Self>>, out: SubmatrixMut<V3, El<Self>>) -> SolveResult
103        where V1: AsPointerToSlice<El<Self>>,
104            V2: AsPointerToSlice<El<Self>>,
105            V3: AsPointerToSlice<El<Self>>
106    {
107        self.get_ring().solve_right(lhs, rhs, out, Global)
108    }
109
110    ///
111    /// Solves a linear system `lhs * X = rhs`.
112    /// 
113    /// For details, see [`LinSolveRing::solve_right()`].
114    /// 
115    #[stability::unstable(feature = "enable")]
116    fn solve_right_with<V1, V2, V3, A>(&self, lhs: SubmatrixMut<V1, El<Self>>, rhs: SubmatrixMut<V2, El<Self>>, out: SubmatrixMut<V3, El<Self>>, allocator: A) -> SolveResult
117        where V1: AsPointerToSlice<El<Self>>,
118            V2: AsPointerToSlice<El<Self>>,
119            V3: AsPointerToSlice<El<Self>>,
120            A: Allocator
121    {
122        self.get_ring().solve_right(lhs, rhs, out, allocator)
123    }
124}
125
126impl<R> LinSolveRingStore for R
127    where R: RingStore, R::Type: LinSolveRing
128{}
129
130impl<R: ?Sized + PrincipalIdealRing> LinSolveRing for R {
131
132    default fn solve_right<V1, V2, V3, A>(&self, lhs: SubmatrixMut<V1, Self::Element>, rhs: SubmatrixMut<V2, Self::Element>, out: SubmatrixMut<V3, Self::Element>, allocator: A) -> SolveResult
133        where V1: AsPointerToSlice<Self::Element>, V2: AsPointerToSlice<Self::Element>, V3: AsPointerToSlice<Self::Element>,
134            A: Allocator
135    {
136        smith::solve_right_using_pre_smith(RingRef::new(self), lhs, rhs, out, allocator)
137    }
138}
139
140impl<R, V, A_ring, C_ring> LinSolveRing for FreeAlgebraImplBase<R, V, A_ring, C_ring>
141    where R: RingStore,
142        R::Type: LinSolveRing,
143        V: VectorView<El<R>>,
144        A_ring: Allocator + Clone, 
145        C_ring: ConvolutionAlgorithm<R::Type>
146{
147    fn solve_right<V1, V2, V3, A>(&self, lhs: SubmatrixMut<V1, Self::Element>, rhs: SubmatrixMut<V2, Self::Element>, out: SubmatrixMut<V3, Self::Element>, allocator: A) -> SolveResult
148        where V1: AsPointerToSlice<Self::Element>,
149            V2: AsPointerToSlice<Self::Element>,
150            V3: AsPointerToSlice<Self::Element>,
151            A:  Allocator
152    {
153        extension::solve_right_over_extension(RingRef::new(self), lhs, rhs, out, allocator)
154    }
155}
feanor_math/algorithms/linsolve/mod.rs

feanor_math/algorithms/linsolve/
mod.rs
