use std::alloc::Global;
use std::cmp::min;
use crate::algorithms::linsolve::smith::determinant_using_pre_smith;
use crate::algorithms::linsolve::LinSolveRing;
use crate::algorithms::linsolve::LinSolveRingStore;
use crate::algorithms::matmul::ComputeInnerProduct;
use crate::algorithms::poly_factor::FactorPolyField;
use crate::divisibility::DivisibilityRing;
use crate::matrix::OwnedMatrix;
use crate::field::*;
use crate::pid::PrincipalIdealRing;
use crate::ring::*;
use crate::seq::*;
use crate::homomorphism::*;
use crate::wrapper::RingElementWrapper;
use super::field::AsField;
use super::field::AsFieldBase;
use super::poly::dense_poly::DensePolyRing;
use super::poly::{PolyRingStore, PolyRing};
pub mod extension_impl;
pub mod galois_field;
pub mod number_field;
pub mod conway;
pub trait FreeAlgebra: RingExtension {
type VectorRepresentation<'a>: VectorFn<El<Self::BaseRing>>
where Self: 'a;
fn canonical_gen(&self) -> Self::Element;
fn rank(&self) -> usize;
fn wrt_canonical_basis<'a>(&'a self, el: &'a Self::Element) -> Self::VectorRepresentation<'a>;
fn from_canonical_basis<V>(&self, vec: V) -> Self::Element
where V: IntoIterator<Item = El<Self::BaseRing>>,
V::IntoIter: DoubleEndedIterator
{
let mut given_len = 0;
let x = self.canonical_gen();
let mut result = self.zero();
for c in vec.into_iter().rev() {
self.mul_assign_ref(&mut result, &x);
self.add_assign(&mut result, self.from(c));
given_len += 1;
}
assert_eq!(given_len, self.rank());
return result;
}
fn from_canonical_basis_extended<V>(&self, vec: V) -> Self::Element
where V: IntoIterator<Item = El<Self::BaseRing>>
{
default_implementations::from_canonical_basis_extended(self, vec)
}
fn charpoly<P, H>(&self, el: &Self::Element, poly_ring: P, hom: H) -> El<P>
where P: RingStore,
P::Type: PolyRing,
<<P::Type as RingExtension>::BaseRing as RingStore>::Type: LinSolveRing,
H: Homomorphism<<Self::BaseRing as RingStore>::Type, <<P::Type as RingExtension>::BaseRing as RingStore>::Type>
{
default_implementations::charpoly(self, el, poly_ring, hom)
}
fn trace(&self, el: Self::Element) -> El<Self::BaseRing> {
let mut current = el;
let generator = self.canonical_gen();
return self.base_ring().sum((0..self.rank()).map(|i| {
let result = self.wrt_canonical_basis(¤t).at(i);
self.mul_assign_ref(&mut current, &generator);
return result;
}));
}
fn discriminant(&self) -> El<Self::BaseRing>
where <Self::BaseRing as RingStore>::Type: PrincipalIdealRing
{
default_implementations::discriminant(self)
}
}
pub trait FreeAlgebraStore: RingStore
where Self::Type: FreeAlgebra
{
delegate!{ FreeAlgebra, fn canonical_gen(&self) -> El<Self> }
delegate!{ FreeAlgebra, fn rank(&self) -> usize }
delegate!{ FreeAlgebra, fn trace(&self, el: El<Self>) -> El<<Self::Type as RingExtension>::BaseRing> }
fn wrt_canonical_basis<'a>(&'a self, el: &'a El<Self>) -> <Self::Type as FreeAlgebra>::VectorRepresentation<'a> {
self.get_ring().wrt_canonical_basis(el)
}
fn from_canonical_basis<V>(&self, vec: V) -> El<Self>
where V: IntoIterator<Item = El<<Self::Type as RingExtension>::BaseRing>>,
V::IntoIter: DoubleEndedIterator
{
self.get_ring().from_canonical_basis(vec)
}
fn from_canonical_basis_extended<V>(&self, vec: V) -> El<Self>
where V: IntoIterator<Item = El<<Self::Type as RingExtension>::BaseRing>>
{
self.get_ring().from_canonical_basis_extended(vec)
}
fn generating_poly<P, H>(&self, poly_ring: P, hom: H) -> El<P>
where P: PolyRingStore,
P::Type: PolyRing,
H: Homomorphism<<<Self::Type as RingExtension>::BaseRing as RingStore>::Type, <<P::Type as RingExtension>::BaseRing as RingStore>::Type>
{
assert!(hom.domain().get_ring() == self.base_ring().get_ring());
poly_ring.sub(
poly_ring.from_terms([(poly_ring.base_ring().one(), self.rank())].into_iter()),
self.poly_repr(&poly_ring, &self.pow(self.canonical_gen(), self.rank()), hom)
)
}
fn as_field(self) -> Result<AsField<Self>, Self>
where Self::Type: DivisibilityRing,
<<Self::Type as RingExtension>::BaseRing as RingStore>::Type: Field + FactorPolyField
{
let poly_ring = DensePolyRing::new(self.base_ring(), "X");
if <_ as FactorPolyField>::factor_poly(&poly_ring, &self.generating_poly(&poly_ring, self.base_ring().identity())).0.len() > 1 {
return Err(self);
} else {
return Ok(RingValue::from(AsFieldBase::promise_is_perfect_field(self)));
}
}
fn poly_repr<P, H>(&self, to: P, el: &El<Self>, hom: H) -> El<P>
where P: PolyRingStore,
P::Type: PolyRing,
H: Homomorphism<<<Self::Type as RingExtension>::BaseRing as RingStore>::Type, <<P::Type as RingExtension>::BaseRing as RingStore>::Type>
{
let coeff_vec = self.wrt_canonical_basis(el);
to.from_terms(
(0..self.rank()).map(|i| coeff_vec.at(i)).enumerate()
.filter(|(_, x)| !self.base_ring().is_zero(x))
.map(|(j, x)| (hom.map(x), j))
)
}
fn discriminant(&self) -> El<<Self::Type as RingExtension>::BaseRing>
where <<Self::Type as RingExtension>::BaseRing as RingStore>::Type: PrincipalIdealRing
{
self.get_ring().discriminant()
}
fn charpoly<P, H>(&self, el: &El<Self>, poly_ring: P, hom: H) -> El<P>
where P: RingStore,
P::Type: PolyRing,
<<P::Type as RingExtension>::BaseRing as RingStore>::Type: LinSolveRing,
H: Homomorphism<<<Self::Type as RingExtension>::BaseRing as RingStore>::Type, <<P::Type as RingExtension>::BaseRing as RingStore>::Type>
{
self.get_ring().charpoly(el, poly_ring, hom)
}
#[stability::unstable(feature = "enable")]
fn with_wrapped_generator<'a, F, const M: usize>(&'a self, f: F) -> [El<Self>; M]
where F: FnOnce(&RingElementWrapper<&'a Self>) -> [RingElementWrapper<&'a Self>; M]
{
let wrapped_indet = RingElementWrapper::new(self, self.canonical_gen());
let mut result_it = f(&wrapped_indet).into_iter();
return std::array::from_fn(|_| result_it.next().unwrap().unwrap());
}
}
mod default_implementations {
use super::*;
pub(super) fn from_canonical_basis_extended<R, V>(ring: &R, vec: V) -> R::Element
where R: ?Sized + FreeAlgebra,
V: IntoIterator<Item = El<<R as RingExtension>::BaseRing>>
{
let mut data = vec.into_iter().collect::<Vec<_>>();
let power_of_canonical_gen = ring.mul(
ring.from_canonical_basis((1..ring.rank()).map(|_| ring.base_ring().zero()).chain([ring.base_ring().one()].into_iter())),
ring.canonical_gen()
);
let mut current_power = ring.one();
return <_ as ComputeInnerProduct>::inner_product(ring, (0..).map_while(|_| {
if data.len() == 0 {
return None;
}
let taken_elements = min(data.len(), ring.rank());
let chunk = data.drain(..taken_elements).chain((taken_elements..ring.rank()).map(|_| ring.base_ring().zero()));
let current = ring.from_canonical_basis(chunk);
let result = (current, ring.clone_el(¤t_power));
ring.mul_assign_ref(&mut current_power, &power_of_canonical_gen);
return Some(result);
}));
}
pub(super) fn charpoly<R, P, H>(ring: &R, el: &R::Element, poly_ring: P, hom: H) -> El<P>
where R: ?Sized + FreeAlgebra,
P: RingStore,
P::Type: PolyRing,
<<P::Type as RingExtension>::BaseRing as RingStore>::Type: LinSolveRing,
H: Homomorphism<<R::BaseRing as RingStore>::Type, <<P::Type as RingExtension>::BaseRing as RingStore>::Type>
{
assert!(!ring.is_zero(el));
let base_ring = hom.codomain();
let mut lhs = OwnedMatrix::zero(ring.rank(), ring.rank(), &base_ring);
let mut current = ring.one();
for j in 0..ring.rank() {
let wrt_basis = ring.wrt_canonical_basis(¤t);
for i in 0..ring.rank() {
*lhs.at_mut(i, j) = hom.map(wrt_basis.at(i));
}
drop(wrt_basis);
ring.mul_assign_ref(&mut current, el);
}
let mut rhs = OwnedMatrix::zero(ring.rank(), 1, &base_ring);
let wrt_basis = ring.wrt_canonical_basis(¤t);
for i in 0..ring.rank() {
*rhs.at_mut(i, 0) = base_ring.negate(hom.map(wrt_basis.at(i)));
}
let mut sol = OwnedMatrix::zero(ring.rank(), 1, &base_ring);
<_ as LinSolveRingStore>::solve_right(base_ring, lhs.data_mut(), rhs.data_mut(), sol.data_mut()).assert_solved();
return poly_ring.from_terms((0..ring.rank()).map(|i| (base_ring.clone_el(sol.at(i, 0)), i)).chain([(base_ring.one(), ring.rank())].into_iter()));
}
pub(super) fn discriminant<R>(ring: &R) -> El<R::BaseRing>
where R: ?Sized + FreeAlgebra,
<R::BaseRing as RingStore>::Type: PrincipalIdealRing
{
let mut current = ring.one();
let generator = ring.canonical_gen();
let traces = (0..(2 * ring.rank())).map(|_| {
let result = ring.trace(ring.clone_el(¤t));
ring.mul_assign_ref(&mut current, &generator);
return result;
}).collect::<Vec<_>>();
let mut matrix = OwnedMatrix::from_fn(ring.rank(), ring.rank(), |i, j| ring.base_ring().clone_el(&traces[i + j]));
let result = determinant_using_pre_smith(ring.base_ring(), matrix.data_mut(), Global);
return result;
}
}
#[stability::unstable(feature = "enable")]
pub fn create_multiplication_matrix<R: FreeAlgebraStore>(ring: R, el: &El<R>) -> OwnedMatrix<El<<R::Type as RingExtension>::BaseRing>>
where R::Type: FreeAlgebra
{
let mut result = OwnedMatrix::zero(ring.rank(), ring.rank(), ring.base_ring());
let mut current = ring.clone_el(el);
let gen = ring.canonical_gen();
for i in 0..ring.rank() {
{
let current_basis_repr = ring.wrt_canonical_basis(¤t);
for j in 0..ring.rank() {
*result.at_mut(j, i) = current_basis_repr.at(j);
}
}
ring.mul_assign_ref(&mut current, &gen);
}
return result;
}
impl<R: RingStore> FreeAlgebraStore for R
where R::Type: FreeAlgebra
{}
#[cfg(any(test, feature = "generic_tests"))]
pub mod generic_tests {
use super::*;
pub fn test_free_algebra_axioms<R: FreeAlgebraStore>(ring: R)
where R::Type: FreeAlgebra
{
let x = ring.canonical_gen();
let n = ring.rank();
let xn_original = ring.pow(ring.clone_el(&x), n);
let xn_vec = ring.wrt_canonical_basis(&xn_original);
let xn = ring.sum(Iterator::map(0..n, |i| ring.mul(ring.inclusion().map(xn_vec.at(i)), ring.pow(ring.clone_el(&x), i))));
assert_el_eq!(ring, xn_original, xn);
let x_n_1_vec_expected = (0..n).map_fn(|i| if i > 0 {
ring.base_ring().add(ring.base_ring().mul(xn_vec.at(n - 1), xn_vec.at(i)), xn_vec.at(i - 1))
} else {
ring.base_ring().mul(xn_vec.at(n - 1), xn_vec.at(0))
});
let x_n_1 = ring.pow(ring.clone_el(&x), n + 1);
let x_n_1_vec_actual = ring.wrt_canonical_basis(&x_n_1);
for i in 0..n {
assert_el_eq!(ring.base_ring(), x_n_1_vec_expected.at(i), x_n_1_vec_actual.at(i));
}
for i in (0..ring.rank()).step_by(5) {
for j in (1..ring.rank()).step_by(7) {
if i == j {
continue;
}
let element = ring.from_canonical_basis((0..n).map(|k| if k == i { ring.base_ring().one() } else if k == j { ring.base_ring().int_hom().map(2) } else { ring.base_ring().zero() }));
let expected = ring.add(ring.pow(ring.clone_el(&x), i), ring.int_hom().mul_map(ring.pow(ring.clone_el(&x), j), 2));
assert_el_eq!(ring, expected, element);
let element_vec = ring.wrt_canonical_basis(&expected);
for k in 0..ring.rank() {
if k == i {
assert_el_eq!(ring.base_ring(), ring.base_ring().one(), element_vec.at(k));
} else if k == j {
assert_el_eq!(ring.base_ring(), ring.base_ring().int_hom().map(2), element_vec.at(k));
} else {
assert_el_eq!(ring.base_ring(), ring.base_ring().zero(), element_vec.at(k));
}
}
}
}
}
}
#[cfg(test)]
use crate::primitive_int::StaticRing;
#[cfg(test)]
use extension_impl::FreeAlgebraImpl;
#[cfg(test)]
use crate::rings::rational::RationalField;
#[test]
fn test_charpoly() {
let ring = FreeAlgebraImpl::new(StaticRing::<i64>::RING, 3, [2, 0, 0]);
let poly_ring = DensePolyRing::new(StaticRing::<i64>::RING, "X");
let [expected] = poly_ring.with_wrapped_indeterminate(|X| [X.pow_ref(3) - 2]);
assert_el_eq!(&poly_ring, &expected, default_implementations::charpoly(ring.get_ring(), &ring.canonical_gen(), &poly_ring, &ring.base_ring().identity()));
let [expected] = poly_ring.with_wrapped_indeterminate(|X| [X.pow_ref(3) - 4]);
assert_el_eq!(&poly_ring, &expected, default_implementations::charpoly(ring.get_ring(), &ring.pow(ring.canonical_gen(), 2), &poly_ring, &ring.base_ring().identity()));
let [expected] = poly_ring.with_wrapped_indeterminate(|X| [X.pow_ref(3) - 6 * X - 6]);
assert_el_eq!(&poly_ring, &expected, default_implementations::charpoly(ring.get_ring(), &ring.add(ring.canonical_gen(), ring.pow(ring.canonical_gen(), 2)), &poly_ring, &ring.base_ring().identity()));
}
#[test]
fn test_trace() {
let ring = FreeAlgebraImpl::new(StaticRing::<i64>::RING, 3, [2, 0, 0]);
assert_eq!(3, ring.trace(ring.from_canonical_basis([1, 0, 0])));
assert_eq!(0, ring.trace(ring.from_canonical_basis([0, 1, 0])));
assert_eq!(0, ring.trace(ring.from_canonical_basis([0, 0, 1])));
assert_eq!(6, ring.trace(ring.from_canonical_basis([2, 0, 0])));
assert_eq!(6, ring.trace(ring.from_canonical_basis([2, 1, 0])));
assert_eq!(6, ring.trace(ring.from_canonical_basis([2, 0, 1])));
}
#[test]
fn test_discriminant() {
let ring = FreeAlgebraImpl::new(StaticRing::<i64>::RING, 3, [2, 0, 0]);
assert_eq!(-108, default_implementations::discriminant(ring.get_ring()));
let ring = FreeAlgebraImpl::new(StaticRing::<i64>::RING, 3, [2, 1, 0]);
assert_eq!(-104, default_implementations::discriminant(ring.get_ring()));
let ring = FreeAlgebraImpl::new(StaticRing::<i64>::RING, 3, [3, 0, 0]);
assert_eq!(-243, default_implementations::discriminant(ring.get_ring()));
let base_ring = DensePolyRing::new(RationalField::new(StaticRing::<i64>::RING), "X");
let [f] = base_ring.with_wrapped_indeterminate(|X| [X.pow_ref(3) + 1]);
let ring = FreeAlgebraImpl::new(&base_ring, 2, [f, base_ring.zero()]);
let [expected] = base_ring.with_wrapped_indeterminate(|X| [4 * X.pow_ref(3) + 4]);
assert_el_eq!(&base_ring, expected, default_implementations::discriminant(ring.get_ring()));
}
#[test]
fn test_from_canonical_basis_extended() {
let ring = FreeAlgebraImpl::new(StaticRing::<i64>::RING, 3, [2]);
let actual = default_implementations::from_canonical_basis_extended(ring.get_ring(), [1, 2, 3, 4, 5, 6, 7]);
let expected = ring.from_canonical_basis([37, 12, 15]);
assert_el_eq!(&ring, expected, actual);
}