feanor-math 3.5.18

use std::alloc::Allocator;
use std::cmp::max;
use std::fmt::{Debug, Formatter, Result};
use std::sync::Arc;

use crate::algorithms;
use crate::algorithms::convolution::ConvolutionAlgorithm;
use crate::algorithms::poly_gcd::PolyTFracGCDRing;
use crate::divisibility::*;
use crate::field::Field;
use crate::integer::{IntegerRing, IntegerRingStore};
use crate::pid::*;
use crate::ring::*;
use crate::rings::poly::*;
use crate::seq::sparse::*;
use crate::seq::{VectorViewMut, *};

/// The univariate polynomial ring `R[X]`. Polynomials are stored as sparse vectors of
/// coefficients, thus giving improved performance in the case that most coefficients are
/// zero.
///
/// Unless polynomials are very sparse, [`crate::rings::poly::dense_poly::DensePolyRing`] will
/// provide better performance.
///
/// # Example
/// ```rust
/// # use feanor_math::ring::*;
/// # use feanor_math::homomorphism::*;
/// # use feanor_math::rings::poly::*;
/// # use feanor_math::rings::poly::sparse_poly::*;
/// # use feanor_math::primitive_int::*;
///
/// let ZZ = StaticRing::<i32>::RING;
/// let P = SparsePolyRing::new(ZZ, "X");
/// let x10_plus_1 = P.add(P.pow(P.indeterminate(), 10), P.int_hom().map(1));
/// let power = P.pow(x10_plus_1, 10);
/// assert_eq!(0, *P.coefficient_at(&power, 1));
/// ```
/// This ring has a [`CanIsoFromTo`] to [`dense_poly::DensePolyRingBase`].
/// ```rust
/// # use feanor_math::assert_el_eq;
/// # use feanor_math::homomorphism::*;
/// # use feanor_math::ring::*;
/// # use feanor_math::rings::poly::*;
/// # use feanor_math::rings::poly::dense_poly::*;
/// # use feanor_math::rings::poly::sparse_poly::*;
/// # use feanor_math::primitive_int::*;
///
/// let ZZ = StaticRing::<i32>::RING;
/// let P = SparsePolyRing::new(ZZ, "X");
/// let P2 = DensePolyRing::new(ZZ, "X");
/// let high_power_of_x = P.pow(P.indeterminate(), 10);
/// assert_el_eq!(
///     P2,
///     P2.pow(P2.indeterminate(), 10),
///     &P.can_iso(&P2).unwrap().map(high_power_of_x)
/// );
/// ```
pub struct SparsePolyRingBase<R: RingStore> {
    base_ring: Arc<R>,
    unknown_name: &'static str,
    zero: El<R>,
}

impl<R: RingStore + Clone> Clone for SparsePolyRingBase<R> {
    fn clone(&self) -> Self {
        SparsePolyRingBase {
            base_ring: self.base_ring.clone(),
            unknown_name: self.unknown_name,
            zero: self.base_ring.zero(),
        }
    }
}

impl<R: RingStore> Debug for SparsePolyRingBase<R>
where
    R::Type: Debug,
{
    fn fmt(&self, f: &mut Formatter<'_>) -> Result {
        f.debug_struct("SparsePolyRingBase")
            .field("base_ring", &self.base_ring.get_ring())
            .finish()
    }
}

/// The univariate polynomial ring `R[X]`, with polynomials being stored as sparse vectors of
/// coefficients. For details, see [`SparsePolyRingBase`].
#[allow(type_alias_bounds)]
pub type SparsePolyRing<R: RingStore> = RingValue<SparsePolyRingBase<R>>;

impl<R: RingStore> SparsePolyRing<R> {
    /// Creates a new sparse polynomial ring over the given base ring.
    pub fn new(base_ring: R, unknown_name: &'static str) -> Self {
        let zero = base_ring.zero();
        Self::from(SparsePolyRingBase {
            base_ring: Arc::new(base_ring),
            unknown_name,
            zero,
        })
    }
}

impl<R: RingStore> SparsePolyRingBase<R> {
    fn degree_truncate(&self, el: &mut SparseMapVector<Arc<R>>) {
        for i in (0..el.len()).rev() {
            if !self.base_ring.is_zero(el.at(i)) {
                el.set_len(i + 1);
                return;
            }
        }
        el.set_len(0);
    }

    fn poly_div<F>(
        &self,
        lhs: &mut <Self as RingBase>::Element,
        rhs: &<Self as RingBase>::Element,
        mut left_div_lc: F,
    ) -> Option<<Self as RingBase>::Element>
    where
        F: FnMut(El<R>) -> Option<El<R>>,
    {
        let lhs_val = std::mem::replace(lhs, self.zero());
        let (quo, rem) = algorithms::poly_div::poly_div_rem(RingRef::new(self), lhs_val, rhs, |x| {
            left_div_lc(self.base_ring().clone_el(x)).ok_or(())
        })
        .ok()?;
        *lhs = rem;
        return Some(quo);
    }
}

/// An element of a [`SparsePolyRing`].
pub struct SparsePolyRingEl<R: RingStore> {
    data: SparseMapVector<Arc<R>>,
}

impl<R: RingStore> Debug for SparsePolyRingEl<R>
where
    El<R>: Debug,
{
    fn fmt(&self, f: &mut Formatter<'_>) -> Result { self.data.fmt(f) }
}

impl<R: RingStore> RingBase for SparsePolyRingBase<R> {
    type Element = SparsePolyRingEl<R>;

    fn clone_el(&self, val: &Self::Element) -> Self::Element { SparsePolyRingEl { data: val.data.clone() } }

    fn add_assign_ref(&self, lhs: &mut Self::Element, rhs: &Self::Element) {
        lhs.data.set_len(max(lhs.data.len(), rhs.data.len()));
        for (i, c) in rhs.data.nontrivial_entries() {
            self.base_ring.add_assign_ref(lhs.data.at_mut(i), c);
        }
        self.degree_truncate(&mut lhs.data);
    }

    fn add_assign(&self, lhs: &mut Self::Element, rhs: Self::Element) { self.add_assign_ref(lhs, &rhs); }

    fn sub_assign_ref(&self, lhs: &mut Self::Element, rhs: &Self::Element) {
        lhs.data.set_len(max(lhs.data.len(), rhs.data.len()));
        for (i, c) in rhs.data.nontrivial_entries() {
            self.base_ring.sub_assign_ref(lhs.data.at_mut(i), c);
        }
        self.degree_truncate(&mut lhs.data);
    }

    fn negate_inplace(&self, lhs: &mut Self::Element) { lhs.data.scan(|_, c| self.base_ring.negate_inplace(c)); }

    fn mul_assign(&self, lhs: &mut Self::Element, rhs: Self::Element) { self.mul_assign_ref(lhs, &rhs); }

    fn mul_assign_ref(&self, lhs: &mut Self::Element, rhs: &Self::Element) { *lhs = self.mul_ref(lhs, rhs); }

    fn zero(&self) -> Self::Element {
        SparsePolyRingEl {
            data: SparseMapVector::new(0, self.base_ring.clone()),
        }
    }

    fn from_int(&self, value: i32) -> Self::Element { self.from(self.base_ring().get_ring().from_int(value)) }

    fn eq_el(&self, lhs: &Self::Element, rhs: &Self::Element) -> bool {
        if lhs.data.len() != rhs.data.len() {
            return false;
        }
        for (i, c) in lhs.data.nontrivial_entries() {
            if !self.base_ring.eq_el(rhs.data.at(i), c) {
                return false;
            }
        }
        for (i, c) in rhs.data.nontrivial_entries() {
            if !self.base_ring.eq_el(lhs.data.at(i), c) {
                return false;
            }
        }
        return true;
    }

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

    fn is_noetherian(&self) -> bool {
        // by Hilbert's basis theorem
        self.base_ring.is_noetherian()
    }

    fn dbg_within<'a>(&self, value: &Self::Element, out: &mut Formatter<'a>, env: EnvBindingStrength) -> Result {
        super::generic_impls::dbg_poly(self, value, out, self.unknown_name, env)
    }

    fn dbg<'a>(&self, value: &Self::Element, out: &mut Formatter<'a>) -> Result {
        self.dbg_within(value, out, EnvBindingStrength::Weakest)
    }

    fn square(&self, value: &mut Self::Element) { *value = self.mul_ref(value, value); }

    fn mul_ref(&self, lhs: &Self::Element, rhs: &Self::Element) -> Self::Element {
        if lhs.data.len() == 0 || rhs.data.len() == 0 {
            return self.zero();
        }
        let mut result = SparseMapVector::new(lhs.data.len() + rhs.data.len() - 1, self.base_ring.clone());
        for (i, c1) in lhs.data.nontrivial_entries() {
            for (j, c2) in rhs.data.nontrivial_entries() {
                self.base_ring
                    .add_assign(result.at_mut(i + j), self.base_ring.mul_ref(c1, c2));
            }
        }
        // if the ring is not zero-divisor free
        self.degree_truncate(&mut result);
        return SparsePolyRingEl { data: result };
    }

    fn mul_assign_int(&self, lhs: &mut Self::Element, rhs: i32) {
        if rhs == 0 {
            *lhs = self.zero();
        } else {
            lhs.data.scan(|_, c| self.base_ring.int_hom().mul_assign_map(c, rhs));
        }
    }

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

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

impl<R> PartialEq for SparsePolyRingBase<R>
where
    R: RingStore,
{
    fn eq(&self, other: &Self) -> bool { self.base_ring.get_ring() == other.base_ring.get_ring() }
}

impl<R: RingStore> RingExtension for SparsePolyRingBase<R> {
    type BaseRing = R;

    fn base_ring(&self) -> &Self::BaseRing { &self.base_ring }

    fn from(&self, x: El<Self::BaseRing>) -> Self::Element {
        let mut result = self.zero();
        if !self.base_ring().is_zero(&x) {
            result.data.set_len(1);
            *result.data.at_mut(0) = x;
        }
        return result;
    }
}

pub trait ImplGenericCanIsoFromToMarker: PolyRing {}

impl<R, A, C> ImplGenericCanIsoFromToMarker for dense_poly::DensePolyRingBase<R, A, C>
where
    R: RingStore,
    A: Allocator + Clone,
    C: ConvolutionAlgorithm<R::Type>,
{
}

impl<R, P> CanHomFrom<P> for SparsePolyRingBase<R>
where
    R: RingStore,
    R::Type: CanHomFrom<<P::BaseRing as RingStore>::Type>,
    P: ImplGenericCanIsoFromToMarker,
{
    type Homomorphism = super::generic_impls::Homomorphism<P, Self>;

    fn has_canonical_hom(&self, from: &P) -> Option<Self::Homomorphism> {
        super::generic_impls::has_canonical_hom(from, self)
    }

    fn map_in(&self, from: &P, el: P::Element, hom: &Self::Homomorphism) -> Self::Element {
        super::generic_impls::map_in(from, self, el, hom)
    }
}

impl<R1, R2> CanHomFrom<SparsePolyRingBase<R1>> for SparsePolyRingBase<R2>
where
    R1: RingStore,
    R2: RingStore,
    R2::Type: CanHomFrom<R1::Type>,
{
    type Homomorphism = <R2::Type as CanHomFrom<R1::Type>>::Homomorphism;

    fn has_canonical_hom(&self, from: &SparsePolyRingBase<R1>) -> Option<Self::Homomorphism> {
        self.base_ring()
            .get_ring()
            .has_canonical_hom(from.base_ring().get_ring())
    }

    fn map_in_ref(
        &self,
        from: &SparsePolyRingBase<R1>,
        el: &SparsePolyRingEl<R1>,
        hom: &Self::Homomorphism,
    ) -> Self::Element {
        let mut result = SparseMapVector::new(el.data.len(), self.base_ring.clone());
        for (j, c) in el.data.nontrivial_entries() {
            *result.at_mut(j) = self
                .base_ring()
                .get_ring()
                .map_in_ref(from.base_ring().get_ring(), c, hom);
        }
        return SparsePolyRingEl { data: result };
    }

    fn map_in(
        &self,
        from: &SparsePolyRingBase<R1>,
        el: <SparsePolyRingBase<R1> as RingBase>::Element,
        hom: &Self::Homomorphism,
    ) -> Self::Element {
        self.map_in_ref(from, &el, hom)
    }
}

impl<R, P> CanIsoFromTo<P> for SparsePolyRingBase<R>
where
    R: RingStore,
    R::Type: CanIsoFromTo<<P::BaseRing as RingStore>::Type>,
    P: ImplGenericCanIsoFromToMarker,
{
    type Isomorphism = super::generic_impls::Isomorphism<P, Self>;

    fn has_canonical_iso(&self, from: &P) -> Option<Self::Isomorphism> {
        self.base_ring()
            .get_ring()
            .has_canonical_iso(from.base_ring().get_ring())
    }

    fn map_out(&self, from: &P, el: Self::Element, iso: &Self::Isomorphism) -> P::Element {
        super::generic_impls::map_out(from, self, el, iso)
    }
}

impl<R1, R2> CanIsoFromTo<SparsePolyRingBase<R1>> for SparsePolyRingBase<R2>
where
    R1: RingStore,
    R2: RingStore,
    R2::Type: CanIsoFromTo<R1::Type>,
{
    type Isomorphism = <R2::Type as CanIsoFromTo<R1::Type>>::Isomorphism;

    fn has_canonical_iso(&self, from: &SparsePolyRingBase<R1>) -> Option<Self::Isomorphism> {
        self.base_ring()
            .get_ring()
            .has_canonical_iso(from.base_ring().get_ring())
    }

    fn map_out(
        &self,
        from: &SparsePolyRingBase<R1>,
        el: Self::Element,
        iso: &Self::Isomorphism,
    ) -> SparsePolyRingEl<R1> {
        let mut result = SparseMapVector::new(el.data.len(), from.base_ring.clone());
        for (j, c) in el.data.nontrivial_entries() {
            *result.at_mut(j) =
                self.base_ring()
                    .get_ring()
                    .map_out(from.base_ring().get_ring(), self.base_ring().clone_el(c), iso);
        }
        return SparsePolyRingEl { data: result };
    }
}

pub struct TermIterator<'a, R>
where
    R: RingStore,
{
    iter: SparseMapVectorIter<'a, Arc<R>>,
}

impl<'a, R> Iterator for TermIterator<'a, R>
where
    R: RingStore,
{
    type Item = (&'a El<R>, usize);

    fn next(&mut self) -> Option<Self::Item> {
        if let Some((i, c)) = self.iter.next() {
            Some((c, i))
        } else {
            None
        }
    }
}

impl<R> PolyRing for SparsePolyRingBase<R>
where
    R: RingStore,
{
    type TermsIterator<'a>
        = TermIterator<'a, R>
    where
        Self: 'a;

    fn indeterminate(&self) -> Self::Element {
        let mut result = self.zero();
        result.data.set_len(2);
        *result.data.at_mut(1) = self.base_ring.one();
        return result;
    }

    fn terms<'a>(&'a self, f: &'a Self::Element) -> TermIterator<'a, R> {
        TermIterator {
            iter: f.data.nontrivial_entries(),
        }
    }

    fn add_assign_from_terms<I>(&self, lhs: &mut Self::Element, rhs: I)
    where
        I: IntoIterator<Item = (El<Self::BaseRing>, usize)>,
    {
        for (c, i) in rhs {
            lhs.data.set_len(max(lhs.data.len(), i + 1));
            self.base_ring().add_assign(lhs.data.at_mut(i), c);
        }
        // if a previously set entry is then set to zero or adds up to zero, this might be not
        // truncated
        self.degree_truncate(&mut lhs.data);
    }

    fn coefficient_at<'a>(&'a self, f: &'a Self::Element, i: usize) -> &'a El<Self::BaseRing> {
        if i < f.data.len() {
            return f.data.at(i);
        } else {
            return &self.zero;
        }
    }

    fn degree(&self, f: &Self::Element) -> Option<usize> { f.data.len().checked_sub(1) }

    fn div_rem_monic(&self, mut lhs: Self::Element, rhs: &Self::Element) -> (Self::Element, Self::Element) {
        assert!(
            self.base_ring()
                .is_one(self.coefficient_at(rhs, self.degree(rhs).unwrap()))
        );
        let quo = self.poly_div(&mut lhs, rhs, Some).unwrap();
        return (quo, lhs);
    }
}

impl<R> Domain for SparsePolyRingBase<R>
where
    R: RingStore,
    R::Type: Domain,
{
}

impl<R> DivisibilityRing for SparsePolyRingBase<R>
where
    R: RingStore,
    R::Type: DivisibilityRing + Domain,
{
    fn checked_left_div(&self, lhs: &Self::Element, rhs: &Self::Element) -> Option<Self::Element> {
        if let Some(d) = self.degree(rhs) {
            let lc = rhs.data.at(d);
            let mut lhs_copy = self.clone_el(lhs);
            let quo = self.poly_div(&mut lhs_copy, rhs, |x| self.base_ring().checked_left_div(&x, lc))?;
            if self.is_zero(&lhs_copy) { Some(quo) } else { None }
        } else if self.is_zero(lhs) {
            Some(self.zero())
        } else {
            None
        }
    }
}

impl<R> PrincipalIdealRing for SparsePolyRingBase<R>
where
    R: RingStore,
    R::Type: Field + PolyTFracGCDRing,
{
    fn checked_div_min(&self, lhs: &Self::Element, rhs: &Self::Element) -> Option<Self::Element> {
        // base ring is a field, so everything is fine
        if self.is_zero(rhs) && self.is_zero(lhs) {
            return Some(self.one());
        } else if self.is_zero(rhs) {
            return None;
        }
        let (quo, rem) = self.euclidean_div_rem(self.clone_el(lhs), rhs);
        if self.is_zero(&rem) {
            return Some(quo);
        } else {
            return None;
        }
    }

    fn extended_ideal_gen(
        &self,
        lhs: &Self::Element,
        rhs: &Self::Element,
    ) -> (Self::Element, Self::Element, Self::Element) {
        algorithms::eea::eea(self.clone_el(lhs), self.clone_el(rhs), RingRef::new(self))
    }

    fn ideal_gen(&self, lhs: &Self::Element, rhs: &Self::Element) -> Self::Element {
        <_ as PolyTFracGCDRing>::gcd(RingRef::new(self), lhs, rhs)
    }
}

impl<R> EuclideanRing for SparsePolyRingBase<R>
where
    R: RingStore,
    R::Type: Field + PolyTFracGCDRing,
{
    fn euclidean_div_rem(&self, mut lhs: Self::Element, rhs: &Self::Element) -> (Self::Element, Self::Element) {
        let lc_inv = self.base_ring.invert(rhs.data.at(self.degree(rhs).unwrap())).unwrap();
        let quo = self
            .poly_div(&mut lhs, rhs, |x| Some(self.base_ring().mul_ref_snd(x, &lc_inv)))
            .unwrap();
        return (quo, lhs);
    }

    fn euclidean_deg(&self, val: &Self::Element) -> Option<usize> {
        return Some(self.degree(val).map(|x| x + 1).unwrap_or(0));
    }
}

#[cfg(test)]
use super::dense_poly::DensePolyRing;
#[cfg(test)]
use crate::primitive_int::StaticRing;
#[cfg(test)]
use crate::rings::finite::FiniteRingStore;
#[cfg(test)]
use crate::rings::zn::zn_static::{Fp, Zn};
#[cfg(test)]
use crate::rings::zn::*;

#[cfg(test)]
fn edge_case_elements<P: PolyRingStore>(poly_ring: P) -> impl Iterator<Item = El<P>>
where
    P::Type: PolyRing,
{
    let base_ring = poly_ring.base_ring();
    vec![
        poly_ring.from_terms([].into_iter()),
        poly_ring.from_terms([(base_ring.int_hom().map(1), 0)].into_iter()),
        poly_ring.from_terms([(base_ring.int_hom().map(1), 1)].into_iter()),
        poly_ring.from_terms([(base_ring.int_hom().map(1), 0), (base_ring.int_hom().map(1), 1)].into_iter()),
        poly_ring.from_terms([(base_ring.int_hom().map(-1), 0)].into_iter()),
        poly_ring.from_terms([(base_ring.int_hom().map(-1), 1)].into_iter()),
        poly_ring.from_terms([(base_ring.int_hom().map(-1), 0), (base_ring.int_hom().map(1), 1)].into_iter()),
        poly_ring.from_terms([(base_ring.int_hom().map(1), 0), (base_ring.int_hom().map(-1), 1)].into_iter()),
        poly_ring.from_terms([(base_ring.int_hom().map(-1), 0), (base_ring.int_hom().map(1), 2)].into_iter()),
        poly_ring.from_terms([(base_ring.int_hom().map(1), 0), (base_ring.int_hom().map(-1), 2)].into_iter()),
        poly_ring.from_terms(
            [
                (base_ring.int_hom().map(1), 0),
                (base_ring.int_hom().map(-1), 2),
                (base_ring.int_hom().map(0), 2),
            ]
            .into_iter(),
        ),
    ]
    .into_iter()
}

#[test]
fn test_ring_axioms() {
    let poly_ring = SparsePolyRing::new(Zn::<7>::RING, "X");
    crate::ring::generic_tests::test_ring_axioms(&poly_ring, edge_case_elements(&poly_ring));
}

#[test]
fn test_poly_ring_axioms() {
    let poly_ring = SparsePolyRing::new(Zn::<7>::RING, "X");
    super::generic_tests::test_poly_ring_axioms(poly_ring, Zn::<7>::RING.elements());
}

#[test]
fn test_canonical_iso_axioms_different_base_ring() {
    let poly_ring1 = SparsePolyRing::new(zn_big::Zn::new(StaticRing::<i128>::RING, 7), "X");
    let poly_ring2 = SparsePolyRing::new(zn_64::Zn::new(7), "X");
    crate::ring::generic_tests::test_hom_axioms(&poly_ring1, &poly_ring2, edge_case_elements(&poly_ring1));
    crate::ring::generic_tests::test_iso_axioms(&poly_ring1, &poly_ring2, edge_case_elements(&poly_ring1));
}

#[test]
fn test_canonical_iso_dense_poly_ring() {
    let poly_ring1 = SparsePolyRing::new(zn_64::Zn::new(7), "X");
    let poly_ring2 = DensePolyRing::new(zn_64::Zn::new(7), "X");
    crate::ring::generic_tests::test_hom_axioms(&poly_ring2, &poly_ring1, edge_case_elements(&poly_ring2));
    crate::ring::generic_tests::test_iso_axioms(&poly_ring2, &poly_ring1, edge_case_elements(&poly_ring2));
}

#[test]
fn test_divisibility_ring_axioms() {
    let poly_ring = SparsePolyRing::new(Fp::<7>::RING, "X");
    crate::divisibility::generic_tests::test_divisibility_axioms(&poly_ring, edge_case_elements(&poly_ring));
}

#[test]
fn test_euclidean_ring_axioms() {
    let poly_ring = SparsePolyRing::new(Fp::<7>::RING, "X");
    crate::pid::generic_tests::test_euclidean_ring_axioms(&poly_ring, edge_case_elements(&poly_ring));
}

#[test]
fn test_principal_ideal_ring_axioms() {
    let poly_ring = SparsePolyRing::new(Fp::<7>::RING, "X");
    crate::pid::generic_tests::test_principal_ideal_ring_axioms(&poly_ring, edge_case_elements(&poly_ring));
}