feanor_math/algorithms/
extension_ops.rs

1use std::alloc::Allocator;
2use std::alloc::Global;
3use std::cmp::min;
4
5use crate::algorithms::int_bisect::bisect_floor;
6use crate::algorithms::linsolve::smith::determinant_using_pre_smith;
7use crate::algorithms::linsolve::*;
8use crate::rings::poly::*;
9use crate::algorithms::matmul::ComputeInnerProduct;
10use crate::divisibility::*;
11use crate::homomorphism::Homomorphism;
12use crate::matrix::OwnedMatrix;
13use crate::pid::PrincipalIdealRing;
14use crate::primitive_int::StaticRing;
15use crate::ring::*;
16use crate::rings::extension::*;
17use crate::seq::*;
18use crate::rings::poly::PolyRing;
19
20///
21/// Default impl for [`FreeAlgebra::from_canonical_basis_extended()`]
22///  
23#[stability::unstable(feature = "enable")]
24pub fn from_canonical_basis_extended<R, V>(ring: &R, vec: V) -> R::Element
25    where R: ?Sized + FreeAlgebra,
26        V: IntoIterator<Item = El<<R as RingExtension>::BaseRing>>
27{
28    let mut data = vec.into_iter().collect::<Vec<_>>();
29    let power_of_canonical_gen = ring.mul(
30        ring.from_canonical_basis((1..ring.rank()).map(|_| ring.base_ring().zero()).chain([ring.base_ring().one()].into_iter())),
31        ring.canonical_gen()
32    );
33    let mut current_power = ring.one();
34    return <_ as ComputeInnerProduct>::inner_product(ring, (0..).map_while(|_| {
35        if data.len() == 0 {
36            return None;
37        }
38        let taken_elements = min(data.len(), ring.rank());
39        let chunk = data.drain(..taken_elements).chain((taken_elements..ring.rank()).map(|_| ring.base_ring().zero()));
40        let current = ring.from_canonical_basis(chunk);
41        let result = (current, ring.clone_el(&current_power));
42        ring.mul_assign_ref(&mut current_power, &power_of_canonical_gen);
43        return Some(result);
44    }));
45}
46
47#[stability::unstable(feature = "enable")]
48pub fn charpoly<R, P, H>(ring: &R, el: &R::Element, poly_ring: P, hom: H) -> El<P>
49    where R: ?Sized + FreeAlgebra,
50        P: RingStore,
51        P::Type: PolyRing,
52        <<P::Type as RingExtension>::BaseRing as RingStore>::Type: LinSolveRing,
53        H: Homomorphism<<R::BaseRing as RingStore>::Type, <<P::Type as RingExtension>::BaseRing as RingStore>::Type>
54{
55    let minpoly = minpoly(ring, el, &poly_ring, hom);
56    let power = StaticRing::<i64>::RING.checked_div(&(ring.rank() as i64), &(poly_ring.degree(&minpoly).unwrap() as i64)).unwrap() as usize;
57    return poly_ring.pow(minpoly, power);
58}
59 
60#[stability::unstable(feature = "enable")]
61pub fn minpoly<R, P, H>(ring: &R, el: &R::Element, poly_ring: P, hom: H) -> El<P>
62    where R: ?Sized + FreeAlgebra,
63        P: RingStore,
64        P::Type: PolyRing,
65        <<P::Type as RingExtension>::BaseRing as RingStore>::Type: LinSolveRing,
66        H: Homomorphism<<R::BaseRing as RingStore>::Type, <<P::Type as RingExtension>::BaseRing as RingStore>::Type>
67{
68    assert!(!ring.is_zero(el));
69    let base_ring = hom.codomain();
70
71    let mut result = None;
72    _ = bisect_floor(StaticRing::<i64>::RING, 1, ring.rank() as i64, |d| {
73        let d = *d as usize;
74        let mut lhs = OwnedMatrix::zero(ring.rank(), d, &base_ring);
75        let mut current = ring.one();
76        for j in 0..d {
77            let wrt_basis = ring.wrt_canonical_basis(&current);
78            for i in 0..ring.rank() {
79                *lhs.at_mut(i, j) = hom.map(wrt_basis.at(i));
80            }
81            drop(wrt_basis);
82            ring.mul_assign_ref(&mut current, el);
83        }
84        let mut rhs = OwnedMatrix::zero(ring.rank(), 1, &base_ring);
85        let wrt_basis = ring.wrt_canonical_basis(&current);
86        for i in 0..ring.rank() {
87            *rhs.at_mut(i, 0) = base_ring.negate(hom.map(wrt_basis.at(i)));
88        }
89        let mut sol = OwnedMatrix::zero(d, 1, &base_ring);
90        let sol_poly = |sol: &OwnedMatrix<El<<P::Type as RingExtension>::BaseRing>>| poly_ring.from_terms((0..d).map(|i| (base_ring.clone_el(sol.at(i, 0)), i)).chain([(base_ring.one(), d)].into_iter()));
91        match <_ as LinSolveRingStore>::solve_right(base_ring, lhs.data_mut(), rhs.data_mut(), sol.data_mut()) {
92            SolveResult::NoSolution => {
93                -1
94            },
95            SolveResult::FoundUniqueSolution => {
96                result = Some(sol_poly(&sol));
97                0
98            },
99            SolveResult::FoundSomeSolution => {
100                result = Some(sol_poly(&sol));
101                1
102            }
103        }
104    });
105
106    return result.unwrap();
107}
108
109///
110/// Default impl for [`FreeAlgebra::discriminant()`]
111/// 
112#[stability::unstable(feature = "enable")]
113pub fn discriminant<R>(ring: &R) -> El<R::BaseRing>
114    where R: ?Sized + FreeAlgebra,
115        <R::BaseRing as RingStore>::Type: PrincipalIdealRing
116{
117    let mut current = ring.one();
118    let generator = ring.canonical_gen();
119    let traces = (0..(2 * ring.rank())).map(|_| {
120        let result = ring.trace(ring.clone_el(&current));
121        ring.mul_assign_ref(&mut current, &generator);
122        return result;
123    }).collect::<Vec<_>>();
124    let mut matrix = OwnedMatrix::from_fn(ring.rank(), ring.rank(), |i, j| ring.base_ring().clone_el(&traces[i + j]));
125    let result = determinant_using_pre_smith(ring.base_ring(), matrix.data_mut(), Global);
126    return result;
127}
128
129#[stability::unstable(feature = "enable")]
130pub fn create_multiplication_matrix<R: RingStore, A: Allocator>(ring: R, el: &El<R>, allocator: A) -> OwnedMatrix<El<<R::Type as RingExtension>::BaseRing>, A>
131    where R::Type: FreeAlgebra
132{
133    let mut result = OwnedMatrix::zero_in(ring.rank(), ring.rank(), ring.base_ring(), allocator);
134    let mut current = ring.clone_el(el);
135    let g = ring.canonical_gen();
136    for i in 0..ring.rank() {
137        {
138            let current_basis_repr = ring.wrt_canonical_basis(&current);
139            for j in 0..ring.rank() {
140                *result.at_mut(j, i) = current_basis_repr.at(j);
141            }
142        }
143        ring.mul_assign_ref(&mut current, &g);
144    }
145    return result;
146}
147
148#[cfg(test)]
149use crate::rings::extension::extension_impl::FreeAlgebraImpl;
150#[cfg(test)]
151use crate::rings::rational::RationalField;
152#[cfg(test)]
153use crate::rings::poly::dense_poly::DensePolyRing;
154
155#[test]
156fn test_charpoly() {
157    let ring = FreeAlgebraImpl::new(StaticRing::<i64>::RING, 3, [2]);
158    let poly_ring = DensePolyRing::new(StaticRing::<i64>::RING, "X");
159
160    let [expected] = poly_ring.with_wrapped_indeterminate(|X| [X.pow_ref(3) - 2]);
161    assert_el_eq!(&poly_ring, &expected, charpoly(ring.get_ring(), &ring.canonical_gen(), &poly_ring, &ring.base_ring().identity()));
162
163    let [expected] = poly_ring.with_wrapped_indeterminate(|X| [X.pow_ref(3) - 4]);
164    assert_el_eq!(&poly_ring, &expected, charpoly(ring.get_ring(), &ring.pow(ring.canonical_gen(), 2), &poly_ring, &ring.base_ring().identity()));
165
166    let [expected] = poly_ring.with_wrapped_indeterminate(|X| [X.pow_ref(3) - 6 * X - 6]);
167    assert_el_eq!(&poly_ring, &expected, charpoly(ring.get_ring(), &ring.add(ring.canonical_gen(), ring.pow(ring.canonical_gen(), 2)), &poly_ring, &ring.base_ring().identity()));
168    
169    let ring = FreeAlgebraImpl::new(StaticRing::<i64>::RING, 4, [2]);
170    let poly_ring = DensePolyRing::new(StaticRing::<i64>::RING, "X");
171
172    let [expected] = poly_ring.with_wrapped_indeterminate(|X| [X.pow_ref(4) - 2]);
173    assert_el_eq!(&poly_ring, &expected, charpoly(ring.get_ring(), &ring.canonical_gen(), &poly_ring, &ring.base_ring().identity()));
174    
175    let [expected] = poly_ring.with_wrapped_indeterminate(|X| [X.pow_ref(4) - 4 * X.pow_ref(2) + 4]);
176    assert_el_eq!(&poly_ring, &expected, charpoly(ring.get_ring(), &ring.pow(ring.canonical_gen(), 2), &poly_ring, &ring.base_ring().identity()));
177}
178
179#[test]
180fn test_minpoly() {
181    let ring = FreeAlgebraImpl::new(StaticRing::<i64>::RING, 6, [2]);
182    let poly_ring = DensePolyRing::new(StaticRing::<i64>::RING, "X");
183
184    let [expected] = poly_ring.with_wrapped_indeterminate(|X| [X.pow_ref(6) - 2]);
185    assert_el_eq!(&poly_ring, &expected, minpoly(ring.get_ring(), &ring.canonical_gen(), &poly_ring, &ring.base_ring().identity()));
186
187    let [expected] = poly_ring.with_wrapped_indeterminate(|X| [X.pow_ref(3) - 2]);
188    assert_el_eq!(&poly_ring, &expected, minpoly(ring.get_ring(), &ring.pow(ring.canonical_gen(), 2), &poly_ring, &ring.base_ring().identity()));
189
190    let [expected] = poly_ring.with_wrapped_indeterminate(|X| [X.pow_ref(2) - 2 * X - 1]);
191    assert_el_eq!(&poly_ring, &expected, minpoly(ring.get_ring(), &ring.add(ring.one(), ring.pow(ring.canonical_gen(), 3)), &poly_ring, &ring.base_ring().identity()));
192}
193
194#[test]
195fn test_trace() {
196    let ring = FreeAlgebraImpl::new(StaticRing::<i64>::RING, 3, [2, 0, 0]);
197
198    assert_eq!(3, ring.trace(ring.from_canonical_basis([1, 0, 0])));
199    assert_eq!(0, ring.trace(ring.from_canonical_basis([0, 1, 0])));
200    assert_eq!(0, ring.trace(ring.from_canonical_basis([0, 0, 1])));
201    assert_eq!(6, ring.trace(ring.from_canonical_basis([2, 0, 0])));
202    assert_eq!(6, ring.trace(ring.from_canonical_basis([2, 1, 0])));
203    assert_eq!(6, ring.trace(ring.from_canonical_basis([2, 0, 1])));
204}
205
206#[test]
207fn test_discriminant() {
208    let ring = FreeAlgebraImpl::new(StaticRing::<i64>::RING, 3, [2, 0, 0]);
209    assert_eq!(-108, discriminant(ring.get_ring()));
210
211    let ring = FreeAlgebraImpl::new(StaticRing::<i64>::RING, 3, [2, 1, 0]);
212    assert_eq!(-104, discriminant(ring.get_ring()));
213    
214    let ring = FreeAlgebraImpl::new(StaticRing::<i64>::RING, 3, [3, 0, 0]);
215    assert_eq!(-243, discriminant(ring.get_ring()));
216
217    let base_ring = DensePolyRing::new(RationalField::new(StaticRing::<i64>::RING), "X");
218    let [f] = base_ring.with_wrapped_indeterminate(|X| [X.pow_ref(3) + 1]);
219    let ring = FreeAlgebraImpl::new(&base_ring, 2, [f, base_ring.zero()]);
220    let [expected] = base_ring.with_wrapped_indeterminate(|X| [4 * X.pow_ref(3) + 4]);
221    assert_el_eq!(&base_ring, expected, discriminant(ring.get_ring()));
222}
223
224#[test]
225fn test_from_canonical_basis_extended() {
226    let ring = FreeAlgebraImpl::new(StaticRing::<i64>::RING, 3, [2]);
227    let actual = from_canonical_basis_extended(ring.get_ring(), [1, 2, 3, 4, 5, 6, 7]);
228    let expected = ring.from_canonical_basis([37, 12, 15]);
229    assert_el_eq!(&ring, expected, actual);
230}
feanor_math/algorithms/extension_ops.rs

feanor_math/algorithms/
extension_ops.rs
