feanor_math/rings/extension/

extension_impl.rs

use std::alloc::Allocator;
use std::alloc::Global;

use serde::de;
use serde::Deserializer;
use serde::Serializer;

use crate::algorithms::convolution::*;
use crate::algorithms::linsolve::LinSolveRing;
use crate::algorithms::poly_factor::FactorPolyField;
use crate::compute_locally::InterpolationBaseRing;
use crate::divisibility::*;
use crate::{impl_localpir_wrap_unwrap_homs, impl_localpir_wrap_unwrap_isos, impl_field_wrap_unwrap_homs, impl_field_wrap_unwrap_isos};
use crate::integer::*;
use crate::iters::multi_cartesian_product;
use crate::iters::MultiProduct;
use crate::matrix::OwnedMatrix;
use crate::primitive_int::StaticRing;
use crate::rings::field::{AsField, AsFieldBase};
use crate::rings::finite::*;
use crate::rings::poly::dense_poly::DensePolyRing;
use crate::seq::*;
use crate::ring::*;
use crate::integer::IntegerRingStore;
use crate::rings::extension::create_multiplication_matrix;
use crate::delegate::DelegateRing;
use crate::homomorphism::*;
use crate::serialization::*;
use crate::seq::sparse::SparseMapVector;

use super::FreeAlgebra;
use super::FreeAlgebraStore;

///
/// Implementation for rings that are generated by a single algebraic element over a base ring,
/// equivalently a quotient of a univariate polynomial ring `R[X]/(f(X))`.
/// 
/// # Example
/// ```
/// # use feanor_math::ring::*;
/// # use feanor_math::rings::extension::*;
/// # use feanor_math::primitive_int::*;
/// # use feanor_math::rings::extension::extension_impl::*;
/// # use feanor_math::assert_el_eq;
/// let base_ring = StaticRing::<i64>::RING;
/// // ring of integers in the 6th cyclotomic number field
/// let ring = FreeAlgebraImpl::new(base_ring, 2, [-1, 1]);
/// let root_of_unity = ring.canonical_gen();
/// assert_el_eq!(ring, ring.one(), ring.pow(root_of_unity, 6));
/// ```
/// 
pub struct FreeAlgebraImplBase<R, V, A = Global, C = KaratsubaAlgorithm>
    where R: RingStore, 
        V: VectorView<El<R>>, 
        A: Allocator + Clone,
        C: ConvolutionAlgorithm<R::Type>
{
    base_ring: R,
    x_pow_rank: V,
    element_allocator: A,
    log2_padded_len: usize,
    rank: usize,
    gen_name: &'static str,
    convolution: C
}

/// 
/// Implementation for rings that are generated by a single algebraic element over a base ring,
/// equivalently a quotient of a univariate polynomial ring `R[X]/(f(X))`. For details, see
/// [`FreeAlgebraImplBase`].
/// 
pub type FreeAlgebraImpl<R, V, A = Global, C = KaratsubaAlgorithm> = RingValue<FreeAlgebraImplBase<R, V, A, C>>;

impl<R, V> FreeAlgebraImpl<R, V>
    where R: RingStore, V: VectorView<El<R>>
{
    ///
    /// Creates a new [`FreeAlgebraImpl`] ring as extension of the given base ring.
    /// 
    /// The created ring is `R[X]/(X^rank - sum_i x_pow_rank[i] X^i)`.
    /// 
    pub fn new(base_ring: R, rank: usize, x_pow_rank: V) -> Self {
        Self::new_with(base_ring, rank, x_pow_rank, "θ", Global, STANDARD_CONVOLUTION)
    }
}

impl<R, V, A, C> Clone for FreeAlgebraImplBase<R, V, A, C>
    where R: RingStore + Clone, 
        V: VectorView<El<R>> + Clone,
        A: Allocator + Clone,
        C: ConvolutionAlgorithm<R::Type> + Clone
{
    fn clone(&self) -> Self {
        Self {
            base_ring: self.base_ring.clone(),
            x_pow_rank: self.x_pow_rank.clone(),
            element_allocator: self.element_allocator.clone(),
            log2_padded_len: self.log2_padded_len,
            rank: self.rank,
            convolution: self.convolution.clone(),
            gen_name: self.gen_name
        }
    }
}

impl<R, V, A, C> Copy for FreeAlgebraImplBase<R, V, A, C>
    where R: RingStore + Copy, 
        V: VectorView<El<R>> + Copy,
        A: Allocator + Copy,
        C: ConvolutionAlgorithm<R::Type> + Copy,
        El<R>: Copy
{}

pub struct FreeAlgebraImplEl<R, A = Global>
    where R: RingStore, A: Allocator
{
    values: Box<[El<R>], A>
}

impl<R, V, A, C> FreeAlgebraImpl<R, V, A, C>
    where R: RingStore, 
        V: VectorView<El<R>>,
        A: Allocator + Clone,
        C: ConvolutionAlgorithm<R::Type>
{
    #[stability::unstable(feature = "enable")]
    pub fn new_with(base_ring: R, rank: usize, x_pow_rank: V, gen_name: &'static str, element_allocator: A, convolution: C) -> Self {
        assert!(rank >= 1);
        assert!(x_pow_rank.len() <= rank);
        let log2_padded_len = StaticRing::<i64>::RING.abs_log2_ceil(&(rank as i64)).unwrap();
        RingValue::from(FreeAlgebraImplBase {
            base_ring, gen_name, x_pow_rank, element_allocator, rank, log2_padded_len, convolution
        })
    }
}

impl<R, V, A, C> FreeAlgebraImplBase<R, V, A, C>
    where R: RingStore, 
        V: VectorView<El<R>>,
        A: Allocator + Clone,
        C: ConvolutionAlgorithm<R::Type>
{
    #[stability::unstable(feature = "enable")]
    pub fn set_convolution<C_new: ConvolutionAlgorithm<R::Type>>(self, new_convolution: C_new) -> FreeAlgebraImpl<R, V, A, C_new> {
        RingValue::from(FreeAlgebraImplBase {
            base_ring: self.base_ring,
            x_pow_rank: self.x_pow_rank,
            element_allocator: self.element_allocator,
            log2_padded_len: self.log2_padded_len,
            rank: self.rank,
            convolution: new_convolution,
            gen_name: self.gen_name
        })
    }

    #[stability::unstable(feature = "enable")]
    pub fn allocator(&self) -> &A {
        &self.element_allocator
    }

    #[stability::unstable(feature = "enable")]
    pub fn gen_name(&self) -> &'static str {
        self.gen_name
    }

    #[stability::unstable(feature = "enable")]
    pub fn x_pow_rank(&self) -> &V {
        &self.x_pow_rank
    }

    #[stability::unstable(feature = "enable")]
    pub fn convolution(&self) -> &C {
        &self.convolution
    }
}

impl<R, V, A, C> PartialEq for FreeAlgebraImplBase<R, V, A, C>
    where R: RingStore, 
        V: VectorView<El<R>>,
        A: Allocator + Clone,
        C: ConvolutionAlgorithm<R::Type>
{
    fn eq(&self, other: &Self) -> bool {
        self.base_ring().get_ring() == other.base_ring().get_ring() && self.rank() == other.rank() &&
            (0..self.rank()).all(|i| (i >= self.x_pow_rank.len() && i >= other.x_pow_rank.len()) ||
                (i >= self.x_pow_rank.len() && self.base_ring().is_zero(other.x_pow_rank.at(i))) ||
                (i >= other.x_pow_rank.len() && self.base_ring().is_zero(self.x_pow_rank.at(i))) ||
                self.base_ring().eq_el(self.x_pow_rank.at(i), other.x_pow_rank.at(i)))
    }
}

impl<R, V, A, C> RingBase for FreeAlgebraImplBase<R, V, A, C>
    where R: RingStore, 
        V: VectorView<El<R>>,
        A: Allocator + Clone,
        C: ConvolutionAlgorithm<R::Type>
{
    type Element = FreeAlgebraImplEl<R, A>;
    
    fn clone_el(&self, val: &Self::Element) -> Self::Element {
        debug_assert_eq!(1 << self.log2_padded_len, val.values.len());
        let mut result_values = Vec::with_capacity_in(val.values.len(), self.element_allocator.clone());
        result_values.extend(val.values.iter().map(|x| self.base_ring().clone_el(x)));
        return FreeAlgebraImplEl {
            values: result_values.into_boxed_slice()
        };
    }

    fn add_assign_ref(&self, lhs: &mut Self::Element, rhs: &Self::Element) {
        debug_assert_eq!(1 << self.log2_padded_len, lhs.values.len());
        debug_assert_eq!(1 << self.log2_padded_len, rhs.values.len());
        for i in 0..self.rank() {
            self.base_ring().add_assign_ref(&mut lhs.values[i], &rhs.values[i]);
        }
    }

    fn add_assign(&self, lhs: &mut Self::Element, rhs: Self::Element) {
        debug_assert_eq!(1 << self.log2_padded_len, lhs.values.len());
        debug_assert_eq!(1 << self.log2_padded_len, rhs.values.len());
        for (i, x) in (0..self.rank()).zip(rhs.values.iter()) {
            self.base_ring().add_assign_ref(&mut lhs.values[i], x);
        }
    }

    fn sub_assign_ref(&self, lhs: &mut Self::Element, rhs: &Self::Element) {
        debug_assert_eq!(1 << self.log2_padded_len, lhs.values.len());
        debug_assert_eq!(1 << self.log2_padded_len, rhs.values.len());
        for i in 0..self.rank() {
            self.base_ring().sub_assign_ref(&mut lhs.values[i], &rhs.values[i]);
        }
    }

    fn negate_inplace(&self, lhs: &mut Self::Element) {
        debug_assert_eq!(1 << self.log2_padded_len, lhs.values.len());
        for i in 0..self.rank() {
            self.base_ring().negate_inplace(&mut lhs.values[i]);
        }
    }

    fn mul_assign(&self, lhs: &mut Self::Element, rhs: Self::Element) {
        debug_assert_eq!(1 << self.log2_padded_len, lhs.values.len());
        debug_assert_eq!(1 << self.log2_padded_len, rhs.values.len());
        self.mul_assign_ref(lhs, &rhs)
    }

    fn mul_assign_ref(&self, lhs: &mut Self::Element, rhs: &Self::Element) {
        debug_assert_eq!(1 << self.log2_padded_len, lhs.values.len());
        debug_assert_eq!(1 << self.log2_padded_len, rhs.values.len());
        let mut tmp = Vec::with_capacity_in(self.rank() * 2, self.element_allocator.clone());
        tmp.extend((0..(2 << self.log2_padded_len)).map(|_| self.base_ring.zero()));
        self.convolution.compute_convolution(&lhs.values[..], &rhs.values[..], &mut tmp[..], self.base_ring());

        struct ReduceModulo<'a, R>
            where R: RingStore
        {
            rank: usize,
            base_ring: &'a R,
            data: &'a mut [El<R>]
        }

        impl<'a, R> SparseVectorViewOperation<El<R>> for ReduceModulo<'a, R>
            where R: RingStore
        {
            type Output<'b> = ()
                where Self: 'b;
        
            fn execute<'b, V: 'b + VectorViewSparse<El<R>>>(self, vector: V) -> () {
                for i in (self.rank..(2 * self.rank)).rev() {
                    for (j, c) in vector.nontrivial_entries() {
                        let add = self.base_ring.mul_ref(c, &mut self.data[i]);
                        self.base_ring.add_assign(&mut self.data[i - self.rank + j], add);
                    }
                }
            }
        }

        let was_sparse = self.x_pow_rank.specialize_sparse(ReduceModulo {
            rank: self.rank,
            base_ring: self.base_ring(),
            data: &mut tmp
        });
        if was_sparse.is_err() {
            for i in (self.rank()..(2 * self.rank())).rev() {
                for j in 0..self.x_pow_rank.len() {
                    let add = self.base_ring.mul_ref(self.x_pow_rank.at(j), &tmp[i]);
                    self.base_ring.add_assign(&mut tmp[i - self.rank() + j], add);
                }
            }
        }
        for i in 0..self.rank() {
            lhs.values[i] = std::mem::replace(&mut tmp[i], self.base_ring.zero());
        }
    }
    
    fn from_int(&self, value: i32) -> Self::Element {
        self.from(self.base_ring().int_hom().map(value))
    }

    fn eq_el(&self, lhs: &Self::Element, rhs: &Self::Element) -> bool {
        debug_assert_eq!(1 << self.log2_padded_len, lhs.values.len());
        debug_assert_eq!(1 << self.log2_padded_len, rhs.values.len());
        (0..self.rank()).all(|i| self.base_ring().eq_el(&lhs.values[i], &rhs.values[i]))
    }

    fn is_zero(&self, value: &Self::Element) -> bool {
        debug_assert_eq!(1 << self.log2_padded_len, value.values.len());
        (0..self.rank()).all(|i| self.base_ring.is_zero(&value.values[i]))
    }

    fn is_one(&self, value: &Self::Element) -> bool {
        debug_assert_eq!(1 << self.log2_padded_len, value.values.len());
        self.base_ring().is_one(&value.values[0]) && (1..self.rank()).all(|i| self.base_ring.is_zero(&value.values[i]))
    }

    fn is_neg_one(&self, value: &Self::Element) -> bool {
        debug_assert_eq!(1 << self.log2_padded_len, value.values.len());
        self.base_ring().is_neg_one(&value.values[0]) && (1..self.rank()).all(|i| self.base_ring.is_zero(&value.values[i]))
    }

    fn is_commutative(&self) -> bool {
        self.base_ring().is_commutative()
    }

    fn is_noetherian(&self) -> bool {
        self.base_ring().is_noetherian()
    }

    fn is_approximate(&self) -> bool {
        self.base_ring().get_ring().is_approximate()
    }

    fn dbg<'a>(&self, value: &Self::Element, out: &mut std::fmt::Formatter<'a>) -> std::fmt::Result {
        debug_assert_eq!(1 << self.log2_padded_len, value.values.len());
        self.dbg_within(value, out, EnvBindingStrength::Weakest)
    }

    fn dbg_within<'a>(&self, value: &Self::Element, out: &mut std::fmt::Formatter<'a>, env: EnvBindingStrength) -> std::fmt::Result {
        debug_assert_eq!(1 << self.log2_padded_len, value.values.len());
        let poly_ring = DensePolyRing::new(self.base_ring(), self.gen_name);
        poly_ring.get_ring().dbg_within(&RingRef::new(self).poly_repr(&poly_ring, value, &self.base_ring().identity()), out, env)
    }

    fn mul_assign_int(&self, lhs: &mut Self::Element, rhs: i32) {
        debug_assert_eq!(1 << self.log2_padded_len, lhs.values.len());
        for i in 0..self.rank() {
            self.base_ring().int_hom().mul_assign_map(&mut lhs.values[i], rhs);
        }
    }

    fn characteristic<I: IntegerRingStore + Copy>(&self, ZZ: I) -> Option<El<I>>
        where I::Type: IntegerRing
    {
        self.base_ring().characteristic(ZZ)
    }
}

impl<R, V, A, C> RingExtension for FreeAlgebraImplBase<R, V, A, C>
    where R: RingStore, 
        V: VectorView<El<R>>,
        A: Allocator + Clone,
        C: ConvolutionAlgorithm<R::Type>
{
    type BaseRing = R;
    
    fn base_ring<'a>(&'a self) -> &'a Self::BaseRing {
        &self.base_ring
    }

    fn from(&self, x: El<Self::BaseRing>) -> Self::Element {
        let mut result_values = Vec::with_capacity_in(1 << self.log2_padded_len, self.element_allocator.clone());
        result_values.push(x);
        result_values.extend((0..((1 << self.log2_padded_len) - 1)).map(|_| self.base_ring().zero()));
        return FreeAlgebraImplEl {
            values: result_values.into_boxed_slice()
        };
    }

    fn mul_assign_base(&self, lhs: &mut Self::Element, rhs: &El<Self::BaseRing>) {
        for i in 0..self.rank() {
            self.base_ring().mul_assign_ref(&mut lhs.values[i], rhs);
        }
    }
}

///
/// It seems like a good idea to also implement [`DivisibilityRing::prepare_divisor()`]
/// and the prepared division functions here. However, currently the interface of
/// [`LinSolveRing`] does not support this. I am not yet sure what is the best approach
/// here, but it seems unlikely that this can be done properly without a breaking release
///  - Either change [`LinSolveRing`] to support "factorizations" of matrices, or rather
///    the equivalent to "prepared divisors" of matrices
///  - Implement [`DivisibilityRing::prepare_divisor()`] only for `R::Type: [`PrincipalIdealRing`]`,
///    which however requires specialization or a specialization workaround
/// 
/// TODO: Decide and implement on the next breaking release
/// 
impl<R, V, A, C> DivisibilityRing for FreeAlgebraImplBase<R, V, A, C>
    where R: RingStore, 
        R::Type: LinSolveRing,
        V: VectorView<El<R>>,
        A: Allocator + Clone,
        C: ConvolutionAlgorithm<R::Type>
{
    fn checked_left_div(&self, lhs: &Self::Element, rhs: &Self::Element) -> Option<Self::Element> {
        let mut mul_matrix: OwnedMatrix<_> = create_multiplication_matrix(RingRef::new(self), rhs);
        let mut lhs_matrix: OwnedMatrix<_> = OwnedMatrix::zero(self.rank(), 1, self.base_ring());
        let data = self.wrt_canonical_basis(&lhs);
        for j in 0..self.rank() {
            *lhs_matrix.at_mut(j, 0) = data.at(j);
        }
        let mut solution: OwnedMatrix<_> = OwnedMatrix::zero(self.rank(), 1, self.base_ring());
        let has_sol = self.base_ring().get_ring().solve_right(mul_matrix.data_mut(), lhs_matrix.data_mut(), solution.data_mut(), Global);
        if has_sol.is_solved() {
            return Some(self.from_canonical_basis((0..self.rank()).map(|i| self.base_ring().clone_el(solution.at(i, 0)))));
        } else {
            return None;
        }
    }

    fn balance_factor<'a, I>(&self, elements: I) -> Option<Self::Element>
        where I: Iterator<Item = &'a Self::Element>, 
            Self: 'a
    {
        self.base_ring().get_ring().balance_factor(elements.flat_map(|x| x.values.iter())).map(|c| RingRef::new(self).inclusion().map(c))
    }

    default fn invert(&self, el: &Self::Element) -> Option<Self::Element> {
        self.checked_left_div(&self.one(), el)
    }
}

impl<R, V, A, C> SerializableElementRing for FreeAlgebraImplBase<R, V, A, C>
    where R: RingStore, 
        V: VectorView<El<R>>,
        A: Allocator + Clone,
        C: ConvolutionAlgorithm<R::Type>,
        R::Type: SerializableElementRing
{
    fn deserialize<'de, D>(&self, deserializer: D) -> Result<Self::Element, D::Error>
        where D: Deserializer<'de>
    {
        let mut data = Vec::with_capacity_in(1 << self.log2_padded_len, self.element_allocator.clone());
        deserialize_seq_helper(deserializer, |x| data.push(x), DeserializeWithRing::new(self.base_ring()))?;
        if data.len() != self.rank() {
            return Err(de::Error::invalid_length(data.len(), &format!("a sequence of {} base ring elements", self.rank()).as_str()));
        }
        data.extend((0..((1 << self.log2_padded_len) - self.rank)).map(|_| self.base_ring().zero()));
        return Ok(FreeAlgebraImplEl { values: data.into_boxed_slice() });
    }

    fn serialize<S>(&self, el: &Self::Element, serializer: S) -> Result<S::Ok, S::Error>
        where S: Serializer
    {
        debug_assert_eq!(1 << self.log2_padded_len, el.values.len());
        serialize_seq_helper(serializer, (&el.values[..self.rank()]).iter().map(|c| SerializeWithRing::new(c, self.base_ring())))
    }
}

impl<R, V, A, C> FreeAlgebra for FreeAlgebraImplBase<R, V, A, C>
    where R: RingStore, 
        V: VectorView<El<R>>,
        A: Allocator + Clone,
        C: ConvolutionAlgorithm<R::Type>
{
    type VectorRepresentation<'a> = CloneElFn<&'a [El<R>], El<R>, CloneRingEl<&'a R>>
        where Self: 'a;

    fn canonical_gen(&self) -> Self::Element {
        if self.rank() == 1 {
            if self.x_pow_rank.len() == 1 {
                self.from_canonical_basis([self.base_ring.clone_el(self.x_pow_rank.at(0))])
            } else {
                assert!(self.x_pow_rank.len() == 0);
                self.zero()
            }
        } else {
            self.from_canonical_basis((0..self.rank()).map(|i| if i == 1 { self.base_ring.one() } else { self.base_ring.zero() }))
        }
    }

    fn wrt_canonical_basis<'a>(&'a self, el: &'a Self::Element) -> Self::VectorRepresentation<'a> {
        (&el.values[..self.rank]).clone_ring_els(self.base_ring())
    }

    fn from_canonical_basis<W>(&self, vec: W) -> Self::Element
        where W: IntoIterator<Item = El<Self::BaseRing>>,
            W::IntoIter: DoubleEndedIterator
    {
        let mut given_len = 0;
        let mut vec_it = vec.into_iter().inspect(|_| given_len += 1);
        let mut result = Vec::with_capacity_in(1 << self.log2_padded_len, self.element_allocator.clone());
        result.extend(vec_it.by_ref());
        result.extend((0..((1 << self.log2_padded_len) - self.rank())).map(|_| self.base_ring().zero()));
        assert!(vec_it.next().is_none());
        assert_eq!(given_len, self.rank());
        FreeAlgebraImplEl { values: result.into_boxed_slice() }
    }

    fn rank(&self) -> usize {
        self.rank
    }
}

impl<R, V, A, C> InterpolationBaseRing for AsFieldBase<FreeAlgebraImpl<R, V, A, C>>
    where R: RingStore + Clone, 
        R::Type: LinSolveRing + FactorPolyField,
        V: VectorView<El<R>>,
        A: Allocator + Clone,
        C: ConvolutionAlgorithm<R::Type>
{
    type ExtendedRingBase<'a> = AsFieldBase<FreeAlgebraImpl<RingRef<'a, Self>, SparseMapVector<RingRef<'a, Self>>, A, KaratsubaAlgorithm<Global>>>
        where Self: 'a;

    type ExtendedRing<'a> = AsField<FreeAlgebraImpl<RingRef<'a, Self>, SparseMapVector<RingRef<'a, Self>>, A, KaratsubaAlgorithm<Global>>>
        where Self: 'a;

    fn in_base<'a, S>(&self, ext_ring: S, el: El<S>) -> Option<Self::Element>
        where Self: 'a, S: RingStore<Type = Self::ExtendedRingBase<'a>>
    {
        let wrt_basis = ext_ring.wrt_canonical_basis(&el);
        if wrt_basis.iter().skip(1).all(|x| self.is_zero(&x)) {
            return Some(wrt_basis.at(0));
        } else {
            return None;
        }
    }

    fn in_extension<'a, S>(&self, ext_ring: S, el: Self::Element) -> El<S>
        where Self: 'a, S: RingStore<Type = Self::ExtendedRingBase<'a>>
    {
        ext_ring.inclusion().map(el)
    }

    fn interpolation_points<'a>(&'a self, count: usize) -> (Self::ExtendedRing<'a>, Vec<El<Self::ExtendedRing<'a>>>) {
        let ZZbig = BigIntRing::RING;
        let characteristic = self.base_ring().characteristic(&ZZbig).unwrap();
        if ZZbig.is_zero(&characteristic) {
            let mut modulus = SparseMapVector::new(1, RingRef::new(self));
            *modulus.at_mut(0) = self.one();
            let field = AsField::from(AsFieldBase::promise_is_perfect_field(FreeAlgebraImpl::new_with(RingRef::new(self), 1, modulus, "X", self.get_delegate().element_allocator.clone(), STANDARD_CONVOLUTION)));
            let points = (0..count).map(|n| field.int_hom().map(n as i32)).collect();
            return (field, points);
        } else {
            unimplemented!()
        }
    }
}

pub struct WRTCanonicalBasisElementCreator<'a, R, V, A, C>
    where R: RingStore, 
        R::Type: FiniteRing, 
        V: VectorView<El<R>>,
        A: Allocator + Clone,
        C: ConvolutionAlgorithm<R::Type>
{
    base_ring: &'a FreeAlgebraImplBase<R, V, A, C>
}

impl<'a, R, V, A, C> Clone for WRTCanonicalBasisElementCreator<'a, R, V, A, C>
    where R: RingStore, 
        R::Type: FiniteRing, 
        V: VectorView<El<R>>,
        A: Allocator + Clone,
        C: ConvolutionAlgorithm<R::Type>
{
    fn clone(&self) -> Self { *self }
}

impl<'a, 'b, R, V, A, C> FnOnce<(&'b [El<R>], )> for WRTCanonicalBasisElementCreator<'a, R, V, A, C>
    where R: RingStore, 
        R::Type: FiniteRing, 
        V: VectorView<El<R>>,
        A: Allocator + Clone,
        C: ConvolutionAlgorithm<R::Type>
{
    type Output = El<FreeAlgebraImpl<R, V, A, C>>;

    extern "rust-call" fn call_once(mut self, args: (&'b [El<R>], )) -> Self::Output {
        self.call_mut(args)
    }
}

impl<'a, 'b, R, V, A, C> FnMut<(&'b [El<R>], )> for WRTCanonicalBasisElementCreator<'a, R, V, A, C>
    where R: RingStore, 
        R::Type: FiniteRing, 
        V: VectorView<El<R>>,
        A: Allocator + Clone,
        C: ConvolutionAlgorithm<R::Type>
{
    extern "rust-call" fn call_mut(&mut self, args: (&'b [El<R>], )) -> Self::Output {
        self.base_ring.from_canonical_basis(args.0.iter().map(|x| self.base_ring.base_ring().clone_el(x)))
    }
}

impl<'a, R, V, A, C> Copy for WRTCanonicalBasisElementCreator<'a, R, V, A, C>
    where R: RingStore, 
        R::Type: FiniteRing, 
        V: VectorView<El<R>>,
        A: Allocator + Clone,
        C: ConvolutionAlgorithm<R::Type>
{}

impl<R, V, A, C> FiniteRing for FreeAlgebraImplBase<R, V, A, C>
    where R: RingStore, 
        R::Type: FiniteRing, 
        V: VectorView<El<R>>,
        A: Allocator + Clone,
        C: ConvolutionAlgorithm<R::Type>
{
    type ElementsIter<'a> = MultiProduct<<R::Type as FiniteRing>::ElementsIter<'a>, WRTCanonicalBasisElementCreator<'a, R, V, A, C>, CloneRingEl<&'a R>, Self::Element>
        where Self: 'a;

    fn elements<'a>(&'a self) -> Self::ElementsIter<'a> {
        multi_cartesian_product((0..self.rank()).map(|_| self.base_ring().elements()), WRTCanonicalBasisElementCreator { base_ring: self }, CloneRingEl(self.base_ring()))
    }

    fn size<I: IntegerRingStore + Copy>(&self, ZZ: I) -> Option<El<I>>
        where I::Type: IntegerRing
    {
        let base_ring_size = self.base_ring().size(ZZ)?;
        if ZZ.get_ring().representable_bits().is_none() || ZZ.abs_log2_ceil(&base_ring_size).unwrap() * self.rank() < ZZ.get_ring().representable_bits().unwrap() {
            Some(ZZ.pow(base_ring_size, self.rank()))
        } else {
            None
        }
    }

    fn random_element<G: FnMut() -> u64>(&self, mut rng: G) -> <Self as RingBase>::Element {
        self.from_canonical_basis((0..self.rank()).map(|_| self.base_ring().random_element(&mut rng)))
    }
}

impl_field_wrap_unwrap_homs!{
    <{R1, V1, A1, C1, R2, V2, A2, C2}> FreeAlgebraImplBase<R1, V1, A1, C1>, FreeAlgebraImplBase<R2, V2, A2, C2>
        where R1: RingStore, R1::Type: LinSolveRing, V1: VectorView<El<R1>>, A1: Allocator + Clone, C1: ConvolutionAlgorithm<R1::Type>,
            R2: RingStore, R2::Type: LinSolveRing, V2: VectorView<El<R2>>, A2: Allocator + Clone, C2: ConvolutionAlgorithm<R2::Type>,
            R2::Type: CanHomFrom<R1::Type>
}

impl_field_wrap_unwrap_isos!{
    <{R1, V1, A1, C1, R2, V2, A2, C2}> FreeAlgebraImplBase<R1, V1, A1, C1>, FreeAlgebraImplBase<R2, V2, A2, C2>
        where R1: RingStore, R1::Type: LinSolveRing, V1: VectorView<El<R1>>, A1: Allocator + Clone, C1: ConvolutionAlgorithm<R1::Type>,
            R2: RingStore, R2::Type: LinSolveRing, V2: VectorView<El<R2>>, A2: Allocator + Clone, C2: ConvolutionAlgorithm<R2::Type>,
            R2::Type: CanIsoFromTo<R1::Type>
}

impl_localpir_wrap_unwrap_homs!{
    <{R1, V1, A1, C1, R2, V2, A2, C2}> FreeAlgebraImplBase<R1, V1, A1, C1>, FreeAlgebraImplBase<R2, V2, A2, C2>
        where R1: RingStore, R1::Type: LinSolveRing, V1: VectorView<El<R1>>, A1: Allocator + Clone, C1: ConvolutionAlgorithm<R1::Type>,
            R2: RingStore, R2::Type: LinSolveRing, V2: VectorView<El<R2>>, A2: Allocator + Clone, C2: ConvolutionAlgorithm<R2::Type>,
            R2::Type: CanHomFrom<R1::Type>
}

impl_localpir_wrap_unwrap_isos!{
    <{R1, V1, A1, C1, R2, V2, A2, C2}> FreeAlgebraImplBase<R1, V1, A1, C1>, FreeAlgebraImplBase<R2, V2, A2, C2>
        where R1: RingStore, R1::Type: LinSolveRing, V1: VectorView<El<R1>>, A1: Allocator + Clone, C1: ConvolutionAlgorithm<R1::Type>,
            R2: RingStore, R2::Type: LinSolveRing, V2: VectorView<El<R2>>, A2: Allocator + Clone, C2: ConvolutionAlgorithm<R2::Type>,
            R2::Type: CanIsoFromTo<R1::Type>
}

impl<R1, V1, A1, C1, R2, V2, A2, C2> CanHomFrom<FreeAlgebraImplBase<R1, V1, A1, C1>> for FreeAlgebraImplBase<R2, V2, A2, C2>
    where R1: RingStore, V1: VectorView<El<R1>>, A1: Allocator + Clone, C1: ConvolutionAlgorithm<R1::Type>,
        R2: RingStore, V2: VectorView<El<R2>>, A2: Allocator + Clone, C2: ConvolutionAlgorithm<R2::Type>,
        R2::Type: CanHomFrom<R1::Type>
{
    type Homomorphism = <R2::Type as CanHomFrom<R1::Type>>::Homomorphism;

    fn has_canonical_hom(&self, from: &FreeAlgebraImplBase<R1, V1, A1, C1>) -> Option<Self::Homomorphism> {
        if self.rank() == from.rank() {
            let hom = self.base_ring.get_ring().has_canonical_hom(from.base_ring.get_ring())?;
            if (0..self.rank()).all(|i| (i >= self.x_pow_rank.len() && i >= from.x_pow_rank.len()) ||
                (i >= self.x_pow_rank.len() && from.base_ring().is_zero(from.x_pow_rank.at(i))) ||
                (i >= from.x_pow_rank.len() && self.base_ring().is_zero(self.x_pow_rank.at(i))) ||
                self.base_ring.eq_el(self.x_pow_rank.at(i), &self.base_ring.get_ring().map_in_ref(from.base_ring.get_ring(), from.x_pow_rank.at(i), &hom))
            ) {
                Some(hom)
            } else {
                None
            }
        } else {
            None
        }
    }

    fn map_in(&self, from: &FreeAlgebraImplBase<R1, V1, A1, C1>, el: <FreeAlgebraImplBase<R1, V1, A1> as RingBase>::Element, hom: &Self::Homomorphism) -> Self::Element {
        self.from_canonical_basis((0..self.rank()).map(|i| self.base_ring.get_ring().map_in_ref(from.base_ring.get_ring(), &el.values[i], hom)))
    }
}

impl<R1, V1, A1, C1, R2, V2, A2, C2> CanIsoFromTo<FreeAlgebraImplBase<R1, V1, A1, C1>> for FreeAlgebraImplBase<R2, V2, A2, C2>
    where R1: RingStore, V1: VectorView<El<R1>>, A1: Allocator + Clone, C1: ConvolutionAlgorithm<R1::Type>,
        R2: RingStore, V2: VectorView<El<R2>>, A2: Allocator + Clone, C2: ConvolutionAlgorithm<R2::Type>,
        R2::Type: CanIsoFromTo<R1::Type>
{
    type Isomorphism = <R2::Type as CanIsoFromTo<R1::Type>>::Isomorphism;

    fn has_canonical_iso(&self, from: &FreeAlgebraImplBase<R1, V1, A1, C1>) -> Option<Self::Isomorphism> {
        if self.rank() == from.rank() {
            let iso = self.base_ring.get_ring().has_canonical_iso(from.base_ring.get_ring())?;
            if (0..self.rank()).all(|i|(i >= self.x_pow_rank.len() && i >= from.x_pow_rank.len()) ||
                (i >= self.x_pow_rank.len() && from.base_ring().is_zero(from.x_pow_rank.at(i))) ||
                (i >= from.x_pow_rank.len() && self.base_ring().is_zero(self.x_pow_rank.at(i))) ||
                from.base_ring.eq_el(&self.base_ring.get_ring().map_out(from.base_ring.get_ring(), self.base_ring.clone_el(self.x_pow_rank.at(i)), &iso), from.x_pow_rank.at(i))
            ) {
                Some(iso)
            } else {
                None
            }
        } else {
            None
        }
    }

    fn map_out(&self, from: &FreeAlgebraImplBase<R1, V1, A1, C1>, el: <FreeAlgebraImplBase<R2, V2, A2> as RingBase>::Element, iso: &Self::Isomorphism) -> <FreeAlgebraImplBase<R1, V1, A1, C1> as RingBase>::Element {
        from.from_canonical_basis((0..self.rank()).map(|i| self.base_ring.get_ring().map_out(from.base_ring.get_ring(), self.base_ring.clone_el(&el.values[i]), iso)))
    }
}

#[cfg(test)]
use crate::rings::zn::zn_64::{Zn, ZnEl};
#[cfg(test)]
use crate::rings::zn::ZnRingStore;
#[cfg(test)]
use crate::rings::zn::zn_static;
#[cfg(test)]
use crate::compute_locally::ToExtRingMap;
#[cfg(test)]
use crate::rings::rational::RationalField;

#[cfg(test)]
fn test_ring0_and_elements() -> (FreeAlgebraImpl<Zn, [ZnEl; 1]>, Vec<FreeAlgebraImplEl<Zn>>) {
    let R = FreeAlgebraImpl::new(Zn::new(7), 1, [Zn::new(7).neg_one()]);
    let elements = R.elements().collect();
    return (R, elements);
}

#[cfg(test)]
fn test_ring1_and_elements() -> (FreeAlgebraImpl<StaticRing::<i64>, [i64; 2]>, Vec<FreeAlgebraImplEl<StaticRing<i64>>>) {
    let ZZ = StaticRing::<i64>::RING;
    let R = FreeAlgebraImpl::new(ZZ, 2, [1, 1]);
    let mut elements = Vec::new();
    for a in -3..=3 {
        for b in -3..=3 {
            elements.push(R.from_canonical_basis([a, b]));
        }
    }
    return (R, elements);
}

#[cfg(test)]
fn test_ring2_and_elements() -> (FreeAlgebraImpl<StaticRing::<i64>, [i64; 3]>, Vec<FreeAlgebraImplEl<StaticRing<i64>>>) {
    let ZZ = StaticRing::<i64>::RING;
    let R = FreeAlgebraImpl::new(ZZ, 3, [1, -1, 1]);
    let mut elements = Vec::new();
    for a in [-2, 0, 1, 2] {
        for b in [-2, 0, 2] {
            for c in [-2, 0, 2] {
                elements.push(R.from_canonical_basis([a, b, c]));
            }
        }
    }
    return (R, elements);
}

#[cfg(test)]
fn test_ring3_and_elements() -> (FreeAlgebraImpl<StaticRing::<i64>, [i64; 1]>, Vec<FreeAlgebraImplEl<StaticRing<i64>>>) {
    let ZZ = StaticRing::<i64>::RING;
    let R = FreeAlgebraImpl::new(ZZ, 3, [1]);
    let mut elements = Vec::new();
    for a in [-2, 0, 1, 2] {
        for b in [-2, 0, 2] {
            for c in [-2, 0, 2] {
                elements.push(R.from_canonical_basis([a, b, c]));
            }
        }
    }
    return (R, elements);
}

#[cfg(test)]
fn test_ring4_and_elements() -> (FreeAlgebraImpl<StaticRing::<i64>, SparseMapVector<StaticRing<i64>>>, Vec<FreeAlgebraImplEl<StaticRing<i64>>>) {
    let ZZ = StaticRing::<i64>::RING;
    // ZZ[X]/(X^3 - X)
    let mut modulus = SparseMapVector::new(3, StaticRing::<i64>::RING);
    *modulus.at_mut(1) = 1;
    let R = FreeAlgebraImpl::new(ZZ, 3, modulus);
    let mut elements = Vec::new();
    for a in [-2, 0, 1, 2] {
        for b in [-2, 0, 2] {
            for c in [-2, 0, 2] {
                elements.push(R.from_canonical_basis([a, b, c]));
            }
        }
    }
    return (R, elements);
}

#[test]
fn test_sparse() {
    let (ring, _) = test_ring4_and_elements();
    assert_el_eq!(ring, ring.canonical_gen(), ring.pow(ring.canonical_gen(), 3));
}

#[test]
fn test_ring_axioms() {
    let (ring, els) = test_ring0_and_elements();
    crate::ring::generic_tests::test_ring_axioms(ring, els.into_iter());
    let (ring, els) = test_ring1_and_elements();
    crate::ring::generic_tests::test_ring_axioms(ring, els.into_iter());
    let (ring, els) = test_ring2_and_elements();
    crate::ring::generic_tests::test_ring_axioms(ring, els.into_iter());
    let (ring, els) = test_ring3_and_elements();
    crate::ring::generic_tests::test_ring_axioms(ring, els.into_iter());
    let (ring, els) = test_ring4_and_elements();
    crate::ring::generic_tests::test_ring_axioms(ring, els.into_iter());

    let base_ring = zn_static::Fp::<257>::RING;
    let x_pow_rank = vec![base_ring.neg_one(); 64];
    let ring = FreeAlgebraImpl::new(base_ring, 64, x_pow_rank);
    let mut rng = oorandom::Rand64::new(0);
    let els = (0..10).map(|_| ring.from_canonical_basis((0..64).map(|_| ring.base_ring().random_element(|| rng.rand_u64()))));
    crate::ring::generic_tests::test_ring_axioms(&ring, els);
}

#[test]
fn test_rank_1_ring() {
    let base_ring = Zn::new(5).as_field().ok().unwrap();
    let ring = FreeAlgebraImpl::new(base_ring, 1, [base_ring.int_hom().map(1)]).as_field().ok().unwrap();
    crate::field::generic_tests::test_field_axioms(&ring, ring.elements());

    assert_el_eq!(ring, ring.one(), ring.canonical_gen());
}

#[test]
fn test_free_algebra_axioms() {
    let (ring, _) = test_ring0_and_elements();
    super::generic_tests::test_free_algebra_axioms(ring);
    let (ring, _) = test_ring1_and_elements();
    super::generic_tests::test_free_algebra_axioms(ring);
    let (ring, _) = test_ring2_and_elements();
    super::generic_tests::test_free_algebra_axioms(ring);
    let (ring, _) = test_ring3_and_elements();
    super::generic_tests::test_free_algebra_axioms(ring);
    let (ring, _) = test_ring4_and_elements();
    super::generic_tests::test_free_algebra_axioms(ring);
}

#[test]
fn test_division() {
    let base_ring = Zn::new(4);
    let i = base_ring.int_hom();
    let ring = FreeAlgebraImpl::new(base_ring, 2, [i.map(-1), i.map(-1)]);

    assert_el_eq!(ring, 
        &ring.from_canonical_basis([i.map(1), i.map(3)]), 
        &ring.checked_div(&ring.one(), &ring.from_canonical_basis([i.map(2), i.map(3)])).unwrap()
    );

    let a = ring.from_canonical_basis([i.map(2), i.map(2)]);
    let b = ring.from_canonical_basis([i.map(0), i.map(2)]);
    assert_el_eq!(ring, a, ring.mul(ring.checked_div(&a, &b).unwrap(), b));

    assert!(ring.checked_div(&ring.one(), &a).is_none());
}

#[test]
fn test_division_ring_of_integers() {
    let base_ring = StaticRing::<i64>::RING;
    let ring = FreeAlgebraImpl::new(base_ring, 2, [11, 0]);

    // the solution to Pell's equation is 10^2 - 3^2 * 11 = 1
    let u = ring.from_canonical_basis([10, 3]);
    let u_inv = ring.from_canonical_basis([10, -3]);

    assert_el_eq!(ring, u_inv, ring.checked_div(&ring.one(), &u).unwrap());
    assert_el_eq!(ring, ring.pow(u_inv, 3), &ring.checked_div(&ring.one(), &ring.pow(u, 3)).unwrap());

    assert!(ring.checked_div(&ring.from_canonical_basis([2, 0]), &ring.from_canonical_basis([3, 0])).is_none());
}

#[test]
fn test_divisibility_axioms() {
    let (ring, els) = test_ring0_and_elements();
    crate::divisibility::generic_tests::test_divisibility_axioms(ring, els.into_iter());
    let (ring, els) = test_ring1_and_elements();
    crate::divisibility::generic_tests::test_divisibility_axioms(ring, els.into_iter());
    let (ring, els) = test_ring2_and_elements();
    crate::divisibility::generic_tests::test_divisibility_axioms(ring, els.into_iter());
    let (ring, els) = test_ring3_and_elements();
    crate::divisibility::generic_tests::test_divisibility_axioms(ring, els.into_iter());
    let (ring, els) = test_ring4_and_elements();
    crate::divisibility::generic_tests::test_divisibility_axioms(ring, els.into_iter());
}

#[test]
fn test_cubic_mul() {
    let base_ring = zn_static::Zn::<256>::RING;
    let modulo = base_ring.int_hom();
    let ring = FreeAlgebraImpl::new(base_ring, 3, [modulo.map(-1), modulo.map(-1), modulo.map(-1)]);
    let a = ring.from_canonical_basis([modulo.map(0), modulo.map(-1), modulo.map(-1)]);
    let b = ring.from_canonical_basis([modulo.map(-1), modulo.map(-2), modulo.map(-1)]);
    assert_el_eq!(ring, b, ring.pow(a, 2));
}

#[test]
fn test_as_field() {
    let base_ring = Zn::new(5).as_field().ok().unwrap();
    let ring = FreeAlgebraImpl::new(base_ring, 1, [base_ring.int_hom().map(1)]).as_field().ok().unwrap();
    crate::field::generic_tests::test_field_axioms(&ring, ring.elements());

    let base_ring = Zn::new(3).as_field().ok().unwrap();
    let ring = FreeAlgebraImpl::new(base_ring, 2, [base_ring.int_hom().map(2)]).as_field().ok().unwrap();
    crate::field::generic_tests::test_field_axioms(&ring, ring.elements());

    assert!(FreeAlgebraImpl::new(base_ring, 2, [base_ring.int_hom().map(1)]).as_field().is_err());
}

#[test]
fn test_serialization() {
    let (ring, els) = test_ring0_and_elements();
    crate::serialization::generic_tests::test_serialization(ring, els.into_iter());
    let (ring, els) = test_ring1_and_elements();
    crate::serialization::generic_tests::test_serialization(ring, els.into_iter());
    let (ring, els) = test_ring2_and_elements();
    crate::serialization::generic_tests::test_serialization(ring, els.into_iter());
    let (ring, els) = test_ring3_and_elements();
    crate::serialization::generic_tests::test_serialization(ring, els.into_iter());
    let (ring, els) = test_ring4_and_elements();
    crate::serialization::generic_tests::test_serialization(ring, els.into_iter());
}

#[test]
#[should_panic]
fn test_from_canonical_basis_enforce_len() {
    let (ring, _) = test_ring1_and_elements();
    _ = ring.from_canonical_basis([0, 1, 2]);
}

#[test]
fn test_interpolation_base_ring() {
    let base_ring = RationalField::new(StaticRing::<i64>::RING);
    let ring = FreeAlgebraImpl::new(base_ring, 3, [base_ring.neg_one(), base_ring.neg_one()]).as_field().ok().unwrap();
    let (ext_map, points) = ToExtRingMap::for_interpolation(ring.get_ring(), 3);
    for _ in points {
        assert_el_eq!(&ring, ring.invert(&ring.canonical_gen()).unwrap(), ext_map.as_base_ring_el(ext_map.codomain().invert(&ext_map.map(ring.canonical_gen())).unwrap()));
    }
}