feanor_math/algorithms/

buchberger.rs

use append_only_vec::AppendOnlyVec;

use crate::computation::*;
use crate::delegate::{UnwrapHom, WrapHom};
use crate::divisibility::{DivisibilityRingStore, DivisibilityRing};
use crate::field::Field;
use crate::homomorphism::Homomorphism;
use crate::local::{PrincipalLocalRing, PrincipalLocalRingStore};
use crate::pid::PrincipalIdealRingStore;
use crate::ring::*;
use crate::seq::*;
use crate::rings::multivariate::*;

use std::cmp::min;
use std::fmt::Debug;

#[stability::unstable(feature = "enable")]
#[derive(PartialEq, Clone, Eq, Hash)]
pub enum SPoly {
    Standard(usize, usize), Nilpotent(/* poly index */ usize, /* power-of-p multiplier */ usize)
}

impl Debug for SPoly {
    fn fmt(&self, f: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
        match self {
            SPoly::Standard(i, j) => write!(f, "S({}, {})", i, j),
            SPoly::Nilpotent(i, k) => write!(f, "p^{} F({})", k, i)
        }
    }
}

fn term_xlcm<P>(ring: P, (l_c, l_m): (&PolyCoeff<P>, &PolyMonomial<P>), (r_c, r_m): (&PolyCoeff<P>, &PolyMonomial<P>)) -> ((PolyCoeff<P>, PolyMonomial<P>), (PolyCoeff<P>, PolyMonomial<P>), (PolyCoeff<P>, PolyMonomial<P>))
    where P: RingStore,
        P::Type: MultivariatePolyRing,
        <<P::Type as RingExtension>::BaseRing as RingStore>::Type: PrincipalLocalRing
{
    let d_c = ring.base_ring().ideal_gen(l_c, r_c);
    let m_m = ring.monomial_lcm(ring.clone_monomial(l_m), r_m);
    let l_factor = ring.base_ring().checked_div(r_c, &d_c).unwrap();
    let r_factor = ring.base_ring().checked_div(l_c, &d_c).unwrap();
    let m_c = ring.base_ring().mul_ref_snd(ring.base_ring().mul_ref_snd(d_c, &r_factor), &l_factor);
    return (
        (l_factor, ring.monomial_div(ring.clone_monomial(&m_m), &l_m).ok().unwrap()),
        (r_factor, ring.monomial_div(ring.clone_monomial(&m_m), &r_m).ok().unwrap()),
        (m_c, m_m)
    );
}

fn term_lcm<P>(ring: P, (l_c, l_m): (&PolyCoeff<P>, &PolyMonomial<P>), (r_c, r_m): (&PolyCoeff<P>, &PolyMonomial<P>)) -> (PolyCoeff<P>, PolyMonomial<P>)
    where P: RingStore,
        P::Type: MultivariatePolyRing,
        <<P::Type as RingExtension>::BaseRing as RingStore>::Type: PrincipalLocalRing
{
    let d_c = ring.base_ring().ideal_gen(l_c, r_c);
    let m_m = ring.monomial_lcm(ring.clone_monomial(l_m), r_m);
    let l_factor = ring.base_ring().checked_div(r_c, &d_c).unwrap();
    let r_factor = ring.base_ring().checked_div(l_c, &d_c).unwrap();
    let m_c = ring.base_ring().mul_ref_snd(ring.base_ring().mul_ref_snd(d_c, &r_factor), &l_factor);
    return (m_c, m_m);
}

impl SPoly {

    #[stability::unstable(feature = "enable")]
    pub fn lcm_term<P, O>(&self, ring: P, basis: &[El<P>], order: O) -> (PolyCoeff<P>, PolyMonomial<P>)
        where P: RingStore + Copy,
            P::Type: MultivariatePolyRing,
            <<P::Type as RingExtension>::BaseRing as RingStore>::Type: PrincipalLocalRing,
            O: MonomialOrder + Copy
    {
        match self {
            SPoly::Standard(i, j) => term_lcm(&ring, ring.LT(&basis[*i], order).unwrap(), ring.LT(&basis[*j], order).unwrap()),
            Self::Nilpotent(i, k) => {
                let (c, m) = ring.LT(&basis[*i], order).unwrap();
                (ring.base_ring().mul_ref_fst(c, ring.base_ring().pow(ring.base_ring().clone_el(ring.base_ring().max_ideal_gen()), *k)), ring.clone_monomial(m))
            }
        }
    }

    
    #[stability::unstable(feature = "enable")]
    pub fn poly<P, O>(&self, ring: P, basis: &[El<P>], order: O) -> El<P>
        where P: RingStore + Copy,
            P::Type: MultivariatePolyRing,
            <<P::Type as RingExtension>::BaseRing as RingStore>::Type: PrincipalLocalRing,
            O: MonomialOrder + Copy
    {
        match self {
            SPoly::Standard(i, j) => {
                let (f1_factor, f2_factor, _) = term_xlcm(&ring, ring.LT(&basis[*i], order).unwrap(), ring.LT(&basis[*j], order).unwrap());
                let mut f1_scaled = ring.clone_el(&basis[*i]);
                ring.mul_assign_monomial(&mut f1_scaled, f1_factor.1);
                ring.inclusion().mul_assign_map(&mut f1_scaled, f1_factor.0);
                let mut f2_scaled = ring.clone_el(&basis[*j]);
                ring.mul_assign_monomial(&mut f2_scaled, f2_factor.1);
                ring.inclusion().mul_assign_map(&mut f2_scaled, f2_factor.0);
                return ring.sub(f1_scaled, f2_scaled);
            },
            SPoly::Nilpotent(i, k) => {
                let mut result = ring.clone_el(&basis[*i]);
                ring.inclusion().mul_assign_map(&mut result, ring.base_ring().pow(ring.base_ring().clone_el(ring.base_ring().max_ideal_gen()), *k));
                return result;
            }
        }
    }
}

#[inline(never)]
fn find_reducer<'a, 'b, P, O, I>(ring: P, f: &El<P>, reducers: I, order: O) -> Option<(usize, &'a El<P>, PolyCoeff<P>, PolyMonomial<P>)>
    where P: RingStore + Copy,
        P::Type: MultivariatePolyRing,
        <<P::Type as RingExtension>::BaseRing as RingStore>::Type: DivisibilityRing,
        O: MonomialOrder + Copy,
        I: Iterator<Item = (&'a El<P>, &'b ExpandedMonomial)>
{
    if ring.is_zero(&f) {
        return None;
    }
    let (f_lc, f_lm) = ring.LT(f, order).unwrap();
    let f_lm_expanded = ring.expand_monomial(f_lm);
    reducers.enumerate().filter_map(|(i, (reducer, reducer_lm_expanded))| {
        if (0..ring.indeterminate_count()).all(|j| reducer_lm_expanded[j] <= f_lm_expanded[j]) {
            let (r_lc, r_lm) = ring.LT(reducer, order).unwrap();
            let quo_m = ring.monomial_div(ring.clone_monomial(f_lm), r_lm).ok().unwrap();
            if let Some(quo_c) = ring.base_ring().checked_div(f_lc, r_lc) {
                return Some((i, reducer, quo_c, quo_m));
            }
        }
        return None;
    }).next()
}

#[inline(never)]
fn filter_spoly<P, O>(ring: P, new_spoly: SPoly, basis: &[El<P>], order: O) -> Option<usize>
    where P: RingStore + Copy,
        P::Type: MultivariatePolyRing,
        <<P::Type as RingExtension>::BaseRing as RingStore>::Type: PrincipalLocalRing,
        O: MonomialOrder + Copy
{
    match new_spoly {
        SPoly::Standard(i, k) => {
            assert!(i < k);
            let (bi_c, bi_m) = ring.LT(&basis[i], order).unwrap();
            let (bk_c, bk_m) = ring.LT(&basis[k], order).unwrap();
            let (S_c, S_m) = term_lcm(ring, (bi_c, bi_m), (bk_c, bk_m));
            let S_c_val = ring.base_ring().valuation(&S_c).unwrap();

            if S_c_val == 0 && order.eq_mon(ring, &ring.monomial_div(ring.clone_monomial(&S_m), &bi_m).ok().unwrap(), &bk_m) {
                return Some(usize::MAX);
            }

            (0..k).filter_map(|j| {
                if j == i {
                    return None;
                }
                // more experiments needed - for some weird reason, replacing "properly divides" with "divides" (assuming
                // I didn't make a mistake) leads to terrible performance
                let (bj_c, bj_m) = ring.LT(&basis[j], order).unwrap();
                let (f_c, f_m) = term_lcm(ring, (bj_c, bj_m), (bk_c, bk_m));
                let f_c_val = ring.base_ring().valuation(&f_c).unwrap();

                if j < i && order.eq_mon(ring, &f_m, &S_m) && f_c_val <= S_c_val {
                    return Some(j);
                }
                if let Ok(quo) = ring.monomial_div(ring.clone_monomial(&S_m), &f_m) {
                    if f_c_val <= S_c_val && (f_c_val < S_c_val || ring.monomial_deg(&quo) > 0) {
                        return Some(j);
                    }
                }
                return None;
            }).next()
        },
        SPoly::Nilpotent(i, k) => {
            let nilpotent_power = ring.base_ring().nilpotent_power().unwrap();
            let f = &basis[i];

            let mut smallest_elim_coeff_valuation = usize::MAX;
            let mut current = ring.LT(f, order).unwrap();
            while ring.base_ring().valuation(&current.0).unwrap() + k >= nilpotent_power {
                smallest_elim_coeff_valuation = min(smallest_elim_coeff_valuation, ring.base_ring().valuation(&current.0).unwrap());
                let next = ring.largest_term_lt(f, order, &current.1);
                if next.is_none() {
                    return Some(usize::MAX);
                }
                current = next.unwrap();
            }
            assert!(smallest_elim_coeff_valuation == usize::MAX || smallest_elim_coeff_valuation + k >= nilpotent_power);
            if smallest_elim_coeff_valuation == usize::MAX || smallest_elim_coeff_valuation + k > nilpotent_power {
                return Some(usize::MAX);
            } else {
                return None;
            }
        }
    }
}

#[stability::unstable(feature = "enable")]
pub fn default_sort_fn<P, O>(ring: P, order: O) -> impl FnMut(&mut [SPoly], &[El<P>])
    where P: RingStore + Copy,
        P::Type: MultivariatePolyRing,
        <<P::Type as RingExtension>::BaseRing as RingStore>::Type: PrincipalLocalRing,
        O: MonomialOrder + Copy
{
    move |open, basis| open.sort_by_key(|spoly| {
        let (lc, lm) = spoly.lcm_term(ring, &basis, order);
        (-(ring.base_ring().valuation(&lc).unwrap_or(0) as i64), -(ring.monomial_deg(&lm) as i64))
    })
}

#[stability::unstable(feature = "enable")]
pub type ExpandedMonomial = Vec<usize>;

fn augment_lm<P, O>(ring: P, f: El<P>, order: O) -> (El<P>, ExpandedMonomial)
    where P: RingStore + Copy,
        P::Type: MultivariatePolyRing,
        <<P::Type as RingExtension>::BaseRing as RingStore>::Type: PrincipalLocalRing,
        O: MonomialOrder
{
    let exponents = ring.expand_monomial(ring.LT(&f, order).unwrap().1);
    return (f, exponents);
}

///
/// Computes a Groebner basis of the ideal generated by the input basis w.r.t. the given term ordering.
/// 
/// For a variant of this function that uses sensible defaults for most parameters, see [`buchberger_simple()`].
/// 
/// The algorithm proceeds F4-style, i.e. reduces multiple S-polynomials before adding them to the basis.
/// When using a fast polynomial ring implementation (e.g. [`crate::rings::multivariate::multivariate_impl::MultivariatePolyRingImpl`]),
/// this makes the algorithm as efficient as standard F4. Furthermore, the behavior can be modified by passing custom functions for
/// `sort_spolys` and `abort_early_if`.
/// 
/// - `sort_spolys` should permute the given list of S-polynomials w.r.t. the given basis; this can be used to customize in which
///   order S-polynomials are reduced, which can have huge impact on performance.
///   Note that S-polynomials that are supposed to be reduced first should be put at the end of the list.
/// - `abort_early_if` takes the current basis (unfortunately, currently with some additional information that can be ignored), and
///   can return `true` to abort the GB computation, yielding the current basis. In this case, the basis will in general not be a GB,
///   but can still be useful (e.g. `abort_early_if` might decide that a GB up to a fixed degree is sufficient).
/// 
/// # Explanation of logging output
/// 
/// If the passed computation controller accepts the logging, it will receive the following symbols:
///  - `-` means an S-polynomial was reduced to zero
///  - `s` means an S-polynomial reduced to a nonzero value and will be added to the basis at the next opportunity
///  - `b(n)` means that the list of all generated basis polynomials has length `n`
///  - `r(n)` means that the current basis of the ideal has length `n`
///  - `S(n)` means that the algorithm still has to reduce `n` more S-polynomials
///  - `f(n)` means that `n` S-polynomials have, in total, been discarded by using the Buchberger criteria
///  - `{n}` means that the algorithm is currently reducing S-polynomials of degree `n`
///  - `!` means that the algorithm decided to discard all current S-polynomial, and restart the computation with the current basis
/// 
#[stability::unstable(feature = "enable")]
pub fn buchberger<P, O, Controller, SortFn, AbortFn>(ring: P, input_basis: Vec<El<P>>, order: O, mut sort_spolys: SortFn, mut abort_early_if: AbortFn, controller: Controller) -> Result<Vec<El<P>>, Controller::Abort>
    where P: RingStore + Copy + Send + Sync,
        El<P>: Send + Sync,
        P::Type: MultivariatePolyRing,
        <P::Type as RingExtension>::BaseRing: Sync,
        <<P::Type as RingExtension>::BaseRing as RingStore>::Type: PrincipalLocalRing,
        O: MonomialOrder + Copy + Send + Sync,
        PolyCoeff<P>: Send + Sync,
        Controller: ComputationController,
        SortFn: FnMut(&mut [SPoly], &[El<P>]),
        AbortFn: FnMut(&[(El<P>, ExpandedMonomial)]) -> bool
{
    controller.run_computation(format_args!("buchberger(len={}, vars={})", input_basis.len(), ring.indeterminate_count()), |controller| {

        // this are the basis polynomials we generated; we only append to this, such that the S-polys remain valid
        let input_basis = inter_reduce(&ring, input_basis.into_iter().map(|f| augment_lm(ring, f, order)).collect(), order).into_iter().map(|(f, _)| f).collect::<Vec<_>>();
        debug_assert!(input_basis.iter().all(|f| !ring.is_zero(f)));

        let nilpotent_power = ring.base_ring().nilpotent_power().and_then(|e| if e != 0 { Some(e) } else { None });
        assert!(nilpotent_power.is_none() || ring.base_ring().is_zero(&ring.base_ring().pow(ring.base_ring().clone_el(ring.base_ring().max_ideal_gen()), nilpotent_power.unwrap())));

        let sort_reducers = |reducers: &mut [(El<P>, ExpandedMonomial)]| {
            // I have no idea why, but this order seems to give the best results
            reducers.sort_by(|(lf, _), (rf, _)| order.compare(ring, &ring.LT(lf, order).unwrap().1, &ring.LT(rf, order).unwrap().1).then_with(|| ring.terms(lf).count().cmp(&ring.terms(rf).count())))
        };

        // invariant: `(reducers) = (basis)` and there exists a reduction to zero for every `f` in `basis` modulo `reducers`;
        // reducers are always stored with an expanded version of their leading monomial, in order to simplify divisibility checks
        let mut reducers: Vec<(El<P>, ExpandedMonomial)> = input_basis.iter().map(|f| augment_lm(ring, ring.clone_el(f), order)).collect::<Vec<_>>();
        sort_reducers(&mut reducers);

        let mut open = Vec::new();
        let mut basis = Vec::new();
        update_basis(ring, input_basis.into_iter(), &mut basis, &mut open, order, nilpotent_power, &mut 0, &mut sort_spolys);

        let mut current_deg = 0;
        let mut filtered_spolys = 0;
        let mut changed = false;
        loop {

            // reduce all known S-polys of minimal lcm degree; in effect, this is the same as 
            // the matrix reduction step during F4
            let spolys_to_reduce_index = open.iter().enumerate().rev().filter(|(_, spoly)| ring.monomial_deg(&spoly.lcm_term(ring, &basis, order).1) > current_deg).next().map(|(i, _)| i + 1).unwrap_or(0);
            let spolys_to_reduce = &open[spolys_to_reduce_index..];

            let computation = ShortCircuitingComputation::new();
            let new_polys = AppendOnlyVec::new();
            let new_polys_ref = &new_polys;
            let basis_ref = &basis;
            let reducers_ref = &reducers;

            computation.handle(controller.clone()).join_many(spolys_to_reduce.as_fn().map_fn(move |spoly| move |handle: ShortCircuitingComputationHandle<(), _>| {
                let mut f = spoly.poly(ring, basis_ref, order);
                
                reduce_poly(ring, &mut f, || reducers_ref.iter().chain(new_polys_ref.iter()).map(|(f, lmf)| (f, lmf)), order);

                if !ring.is_zero(&f) {
                    log_progress!(handle, "s");
                    _ = new_polys_ref.push(augment_lm(ring, f, order));
                } else {
                    log_progress!(handle, "-");
                }

                checkpoint!(handle);
                return Ok(None);
            }));

            drop(open.drain(spolys_to_reduce_index..));
            let new_polys = new_polys.into_vec();
            _ = computation.finish()?;

            // process the generated new polynomials
            if new_polys.len() == 0 && open.len() == 0 {
                if changed {
                    log_progress!(controller, "!");
                    // this seems necessary, as the invariants for `reducers` don't imply that it already is a GB;
                    // more concretely, reducers contains polys of basis that are reduced with eath other, but the
                    // S-polys between two of them might not have been considered
                    return buchberger::<P, O, _, _, _>(ring, reducers.into_iter().map(|(f, _)| f).collect(), order, sort_spolys, abort_early_if, controller);
                } else {
                    return Ok(reducers.into_iter().map(|(f, _)| f).collect());
                }
            } else if new_polys.len() == 0 {
                current_deg = ring.monomial_deg(&open.last().unwrap().lcm_term(ring, &basis, order).1);
                log_progress!(controller, "{{{}}}", current_deg);
            } else {
                changed = true;
                current_deg = 0;
                update_basis(ring, new_polys.iter().map(|(f, _)| ring.clone_el(f)), &mut basis, &mut open, order, nilpotent_power, &mut filtered_spolys, &mut sort_spolys);
                log_progress!(controller, "(b={})(S={})(f={})", basis.len(), open.len(), filtered_spolys);

                reducers.extend(new_polys.into_iter());
                reducers = inter_reduce(ring, reducers, order);
                sort_reducers(&mut reducers);
                log_progress!(controller, "(r={})", reducers.len());
                if abort_early_if(&reducers) {
                    log_progress!(controller, "(early_abort)");
                    return Ok(reducers.into_iter().map(|(f, _)| f).collect());
                }
            }

            // less S-polys if we restart from scratch with reducers
            if open.len() + filtered_spolys > reducers.len() * reducers.len() / 2 + reducers.len() * nilpotent_power.unwrap_or(0) + 1 {
                log_progress!(controller, "!");
                return buchberger::<P, O, _, _, _>(ring, reducers.into_iter().map(|(f, _)| f).collect(), order, sort_spolys, abort_early_if, controller);
            }
        }
    })
}

fn update_basis<I, P, O, SortFn>(ring: P, new_polys: I, basis: &mut Vec<El<P>>, open: &mut Vec<SPoly>, order: O, nilpotent_power: Option<usize>, filtered_spolys: &mut usize, sort_spolys: &mut SortFn)
    where P: RingStore + Copy,
        P::Type: MultivariatePolyRing,
        <<P::Type as RingExtension>::BaseRing as RingStore>::Type: PrincipalLocalRing,
        O: MonomialOrder + Copy,
        SortFn: FnMut(&mut [SPoly], &[El<P>]),
        I: Iterator<Item = El<P>>
{
    for new_poly in new_polys {
        basis.push(new_poly);
        for i in 0..(basis.len() - 1) {
            let spoly = SPoly::Standard(i, basis.len() - 1);
            if filter_spoly(ring, spoly.clone(), &*basis, order).is_none() {
                open.push(spoly);
            } else {
                *filtered_spolys += 1;
            }
        }
        if let Some(e) = nilpotent_power {
            for k in 1..e {
                let spoly = SPoly::Nilpotent(basis.len() - 1, k);
                if filter_spoly(ring, spoly.clone(), &*basis, order).is_none() {
                    open.push(spoly);
                } else {
                    *filtered_spolys += 1;
                }
            }
        }
    }
    sort_spolys(&mut *open, &*basis);
}

fn reduce_poly<'a, 'b, F, I, P, O>(ring: P, to_reduce: &mut El<P>, mut reducers: F, order: O)
    where P: 'a + RingStore + Copy,
        P::Type: MultivariatePolyRing,
        <<P::Type as RingExtension>::BaseRing as RingStore>::Type: DivisibilityRing,
        O: MonomialOrder + Copy,
        F: FnMut() -> I,
        I: Iterator<Item = (&'a El<P>, &'b ExpandedMonomial)>
{
    while let Some((_, reducer, quo_c, quo_m)) = find_reducer(ring, to_reduce, reducers(), order) {
        let prev_lm = ring.clone_monomial(&ring.LT(to_reduce, order).unwrap().1);
        let mut scaled_reducer = ring.clone_el(reducer);
        ring.mul_assign_monomial(&mut scaled_reducer, ring.clone_monomial(&quo_m));
        ring.inclusion().mul_assign_ref_map(&mut scaled_reducer, &quo_c);
        debug_assert!(order.compare(ring, ring.LT(&scaled_reducer, order).unwrap().1, ring.LT(&to_reduce, order).unwrap().1) == std::cmp::Ordering::Equal);
        ring.sub_assign(to_reduce, scaled_reducer);
        debug_assert!(ring.is_zero(&to_reduce) || order.compare(ring, &ring.LT(&to_reduce, order).unwrap().1, &prev_lm) == std::cmp::Ordering::Less);
    }
}

#[stability::unstable(feature = "enable")]
pub fn multivariate_division<'a, P, O, I>(ring: P, mut f: El<P>, reducers: I, order: O) -> El<P>
    where P: 'a + RingStore + Copy,
        P::Type: MultivariatePolyRing,
        <<P::Type as RingExtension>::BaseRing as RingStore>::Type: DivisibilityRing,
        O: MonomialOrder + Copy,
        I: Clone + Iterator<Item = &'a El<P>>
{
    let lms = reducers.clone().map(|f| ring.expand_monomial(ring.LT(f, order).unwrap().1)).collect::<Vec<_>>();
    reduce_poly(ring, &mut f, || reducers.clone().zip(lms.iter()), order);
    return f;
}

fn inter_reduce<P, O>(ring: P, mut polys: Vec<(El<P>, ExpandedMonomial)>, order: O) -> Vec<(El<P>, ExpandedMonomial)>
    where P: RingStore + Copy,
        P::Type: MultivariatePolyRing,
        <<P::Type as RingExtension>::BaseRing as RingStore>::Type: DivisibilityRing,
        O: MonomialOrder + Copy
{
    let mut changed = true;
    while changed {
        changed = false;
        let mut i = 0;
        while i < polys.len() {
            let last_i = polys.len() - 1;
            polys.swap(i, last_i);
            let (reducers, to_reduce) = polys.split_at_mut(last_i);
            let to_reduce = &mut to_reduce[0];

            reduce_poly(ring, &mut to_reduce.0, || reducers.iter().map(|(f, lmf)| (f, lmf)), order);

            // undo swap so that the outer loop still iterates over every poly
            if !ring.is_zero(&to_reduce.0) {
                to_reduce.1 = ring.expand_monomial(ring.LT(&to_reduce.0, order).unwrap().1);
                polys.swap(i, last_i);
                i += 1;
            } else {
                _ = polys.pop().unwrap();
            }
        }
    }
    return polys;
}

use crate::rings::local::AsLocalPIR;
use crate::rings::multivariate::multivariate_impl::MultivariatePolyRingImpl;

///
/// Computes a Groebner basis of the ideal generated by the input basis w.r.t. the given term ordering.
/// 
/// For a variant of this function that allows for more configuration, see [`buchberger()`].
/// 
pub fn buchberger_simple<P, O>(ring: P, input_basis: Vec<El<P>>, order: O) -> Vec<El<P>>
    where P: RingStore + Copy + Send + Sync,
        El<P>: Send + Sync,
        P::Type: MultivariatePolyRing,
        <P::Type as RingExtension>::BaseRing: Sync,
        <<P::Type as RingExtension>::BaseRing as RingStore>::Type: Field,
        O: MonomialOrder + Copy + Send + Sync,
        PolyCoeff<P>: Send + Sync
{
    let as_local_pir = AsLocalPIR::from_field(ring.base_ring());
    let new_poly_ring = MultivariatePolyRingImpl::new(&as_local_pir, ring.indeterminate_count());
    let from_ring = new_poly_ring.lifted_hom(ring, WrapHom::to_delegate_ring(as_local_pir.get_ring()));
    let result = buchberger::<_, _, _, _, _>(
        &new_poly_ring, 
        input_basis.into_iter().map(|f| from_ring.map(f)).collect(), 
        order, 
        default_sort_fn(&new_poly_ring, order), 
        |_| false, 
        DontObserve
    ).unwrap_or_else(no_error);
    let to_ring = ring.lifted_hom(&new_poly_ring, UnwrapHom::from_delegate_ring(as_local_pir.get_ring()));
    return result.into_iter().map(|f| to_ring.map(f)).collect();
}

#[cfg(test)]
use crate::rings::poly::{dense_poly, PolyRingStore};
#[cfg(test)]
use crate::rings::zn::zn_static;
#[cfg(test)]
use crate::integer::BigIntRing;
#[cfg(test)]
use crate::rings::rational::RationalField;

#[test]
fn test_buchberger_small() {
    let base = zn_static::F17;
    let ring = MultivariatePolyRingImpl::new(base, 2);

    let f1 = ring.from_terms([
        (1, ring.create_monomial([2, 0])),
        (1, ring.create_monomial([0, 2])),
        (16, ring.create_monomial([0, 0]))
    ].into_iter());
    let f2 = ring.from_terms([
        (1, ring.create_monomial([1, 1])),
        (15, ring.create_monomial([0, 0]))
    ].into_iter());

    let actual = buchberger(&ring, vec![ring.clone_el(&f1), ring.clone_el(&f2)], DegRevLex, default_sort_fn(&ring, DegRevLex), |_| false, TEST_LOG_PROGRESS).unwrap_or_else(no_error);

    let expected = ring.from_terms([
        (16, ring.create_monomial([0, 3])),
        (15, ring.create_monomial([1, 0])),
        (1, ring.create_monomial([0, 1])),
    ].into_iter());

    assert_eq!(3, actual.len());
    assert_el_eq!(ring, ring.zero(), multivariate_division(&ring, f1, actual.iter(), DegRevLex));
    assert_el_eq!(ring, ring.zero(), multivariate_division(&ring, f2, actual.iter(), DegRevLex));
    assert_el_eq!(ring, ring.zero(), multivariate_division(&ring, expected, actual.iter(), DegRevLex));
}

#[test]
fn test_buchberger_larger() {
    let base = zn_static::F17;
    let ring = MultivariatePolyRingImpl::new(base, 3);

    let f1 = ring.from_terms([
        (1, ring.create_monomial([2, 1, 1])),
        (1, ring.create_monomial([0, 2, 0])),
        (1, ring.create_monomial([1, 0, 1])),
        (2, ring.create_monomial([1, 0, 0])),
        (1, ring.create_monomial([0, 0, 0]))
    ].into_iter());
    let f2 = ring.from_terms([
        (1, ring.create_monomial([0, 3, 1])),
        (1, ring.create_monomial([0, 0, 3])),
        (1, ring.create_monomial([1, 1, 0]))
    ].into_iter());
    let f3 = ring.from_terms([
        (1, ring.create_monomial([1, 0, 2])),
        (1, ring.create_monomial([1, 0, 1])),
        (2, ring.create_monomial([0, 1, 1])),
        (7, ring.create_monomial([0, 0, 0]))
    ].into_iter());

    let actual = buchberger(&ring, vec![ring.clone_el(&f1), ring.clone_el(&f2), ring.clone_el(&f3)], DegRevLex, default_sort_fn(&ring, DegRevLex), |_| false, TEST_LOG_PROGRESS).unwrap_or_else(no_error);

    let g1 = ring.from_terms([
        (1, ring.create_monomial([0, 4, 0])),
        (8, ring.create_monomial([0, 3, 1])),
        (12, ring.create_monomial([0, 1, 3])),
        (6, ring.create_monomial([0, 0, 4])),
        (1, ring.create_monomial([0, 3, 0])),
        (13, ring.create_monomial([0, 2, 1])),
        (11, ring.create_monomial([0, 1, 2])),
        (10, ring.create_monomial([0, 0, 3])),
        (11, ring.create_monomial([0, 2, 0])),
        (12, ring.create_monomial([0, 1, 1])),
        (6, ring.create_monomial([0, 0, 2])),
        (6, ring.create_monomial([0, 1, 0])),
        (13, ring.create_monomial([0, 0, 1])),
        (9, ring.create_monomial([0, 0, 0]))
    ].into_iter());

    assert_el_eq!(ring, ring.zero(), multivariate_division(&ring, g1, actual.iter(), DegRevLex));
}

#[test]
fn test_generic_computation() {
    let base = zn_static::F17;
    let ring = MultivariatePolyRingImpl::new(base, 6);
    let poly_ring = dense_poly::DensePolyRing::new(&ring, "X");

    let var_i = |i: usize| ring.create_term(base.one(), ring.create_monomial((0..ring.indeterminate_count()).map(|j| if i == j { 1 } else { 0 })));
    let X1 = poly_ring.mul(
        poly_ring.from_terms([(var_i(0), 0), (ring.one(), 1)].into_iter()),
        poly_ring.from_terms([(var_i(1), 0), (ring.one(), 1)].into_iter())
    );
    let X2 = poly_ring.mul(
        poly_ring.add(poly_ring.clone_el(&X1), poly_ring.from_terms([(var_i(2), 0), (var_i(3), 1)].into_iter())),
        poly_ring.add(poly_ring.clone_el(&X1), poly_ring.from_terms([(var_i(4), 0), (var_i(5), 1)].into_iter()))
    );
    let basis = vec![
        ring.sub_ref_snd(ring.int_hom().map(1), poly_ring.coefficient_at(&X2, 0)),
        ring.sub_ref_snd(ring.int_hom().map(1), poly_ring.coefficient_at(&X2, 1)),
        ring.sub_ref_snd(ring.int_hom().map(1), poly_ring.coefficient_at(&X2, 2)),
    ];

    let start = std::time::Instant::now();
    let gb1 = buchberger(&ring, basis.iter().map(|f| ring.clone_el(f)).collect(), DegRevLex, default_sort_fn(&ring, DegRevLex), |_| false, TEST_LOG_PROGRESS).unwrap_or_else(no_error);
    let end = std::time::Instant::now();

    println!("Computed GB in {} ms", (end - start).as_millis());

    assert_eq!(11, gb1.len());
}

#[test]
fn test_gb_local_ring() {
    let base = AsLocalPIR::from_zn(zn_static::Zn::<16>::RING).unwrap();
    let ring: MultivariatePolyRingImpl<_> = MultivariatePolyRingImpl::new(base, 1);
    
    let f = ring.from_terms([(base.int_hom().map(4), ring.create_monomial([1])), (base.one(), ring.create_monomial([0]))].into_iter());
    let gb = buchberger(&ring, vec![f], DegRevLex, default_sort_fn(&ring, DegRevLex), |_| false, TEST_LOG_PROGRESS).unwrap_or_else(no_error);

    assert_eq!(1, gb.len());
    assert_el_eq!(ring, ring.one(), gb[0]);
}

#[test]
fn test_gb_lex() {
    let ZZ = BigIntRing::RING;
    let QQ = AsLocalPIR::from_field(RationalField::new(ZZ));
    let QQYX = MultivariatePolyRingImpl::new(&QQ, 2);
    let [f, g] = QQYX.with_wrapped_indeterminates(|[Y, X]| [ 1 + X.pow_ref(2) + 2 * Y + (1 + X) * Y.pow_ref(2), 3 + X + (2 + X) * Y + (1 + X + X.pow_ref(2)) * Y.pow_ref(2) ]);
    let expected = QQYX.with_wrapped_indeterminates(|[Y, X]| [ 
        X.pow_ref(8) + 2 * X.pow_ref(7) + 3 * X.pow_ref(6) - 5 * X.pow_ref(5) - 10 * X.pow_ref(4) - 7 * X.pow_ref(3) + 8 * X.pow_ref(2) + 8 * X + 4,
        2 * Y + X.pow_ref(6) + 3 * X.pow_ref(5) + 6 * X.pow_ref(4) + X.pow_ref(3) - 7 * X.pow_ref(2) - 12 * X - 2, 
    ]);

    let mut gb = buchberger_simple(&QQYX, vec![f, g], Lex);

    assert_eq!(2, gb.len());
    gb.sort_unstable_by_key(|f| QQYX.appearing_indeterminates(f).len());
    for (mut f, mut e) in gb.into_iter().zip(expected.into_iter()) {
        let f_lc_inv = QQ.invert(QQYX.LT(&f, Lex).unwrap().0).unwrap();
        QQYX.inclusion().mul_assign_map(&mut f, f_lc_inv);
        let e_lc_inv = QQ.invert(QQYX.LT(&e, Lex).unwrap().0).unwrap();
        QQYX.inclusion().mul_assign_map(&mut e, e_lc_inv);
        assert_el_eq!(QQYX, e, f);
    }
}

#[cfg(test)]
#[cfg(feature = "parallel")]
static TEST_COMPUTATION_CONTROLLER: ExecuteMultithreaded<LogProgress> = RunMultithreadedLogProgress;
#[cfg(test)]
#[cfg(not(feature = "parallel"))]
static TEST_COMPUTATION_CONTROLLER: LogProgress = TEST_LOG_PROGRESS;

#[ignore]
#[test]
fn test_expensive_gb_1() {
    let base = AsLocalPIR::from_zn(zn_static::Zn::<16>::RING).unwrap();
    let ring: MultivariatePolyRingImpl<_> = MultivariatePolyRingImpl::new(base, 12);

    let system = ring.with_wrapped_indeterminates(|[Y0, Y1, Y2, Y3, Y4, Y5, Y6, Y7, Y8, Y9, Y10, Y11]| [
        Y0 * Y1 * Y2 * Y3.pow_ref(2) * Y4.pow_ref(2) + 4 * Y0 * Y1 * Y2 * Y3 * Y4 * Y4 * Y8 + Y0 * Y1 * Y2 * Y5.pow_ref(2) * Y8.pow_ref(2) + Y0 * Y2 * Y3 * Y4 * Y6 + Y0 * Y1 * Y3 * Y4 * Y7 + Y0 * Y2 * Y5 * Y6 * Y8 + Y0 * Y1 * Y5 * Y7 * Y8 + Y0 * Y2 * Y3 * Y5 * Y10 + Y0 * Y1 * Y3 * Y5 * Y11 + Y0 * Y6 * Y7 + Y3 * Y5 * Y9 -  4,
        2 * Y0 * Y1 * Y2 * Y3.pow_ref(2) * Y4 * Y5 + 2 * Y0 * Y1 * Y2 * Y3 * Y5.pow_ref(2) * Y8 + Y0 * Y2 * Y3 * Y5 * Y6 + Y0 * Y1 * Y3 * Y5 * Y7 + 8,
        Y0 * Y1 * Y2 * Y3.pow_ref(2) * Y5.pow_ref(2) - 5
    ]);

    let part_of_result = ring.with_wrapped_indeterminates(|[_Y0, Y1, Y2, _Y3, _Y4, _Y5, Y6, Y7, _Y8, _Y9, _Y10, _Y11]| [
        4 * Y2.pow_ref(2) * Y6.pow_ref(2) -  4 * Y1.pow_ref(2) * Y7.pow_ref(2),
        8 * Y2 * Y6 + 8 * Y1 * Y7.clone()
    ]);

    let start = std::time::Instant::now();
    let gb = buchberger(&ring, system.iter().map(|f| ring.clone_el(f)).collect(), DegRevLex, default_sort_fn(&ring, DegRevLex), |_| false, TEST_COMPUTATION_CONTROLLER).unwrap_or_else(no_error);
    let end = std::time::Instant::now();

    println!("Computed GB in {} ms", (end - start).as_millis());

    for f in &part_of_result {
        assert!(ring.is_zero(&multivariate_division(&ring, ring.clone_el(f), gb.iter(), DegRevLex)));
    }

    assert_eq!(108, gb.len());
}

#[test]
#[ignore]
fn test_expensive_gb_2() {
    let base = zn_static::Fp::<7>::RING;
    let ring = MultivariatePolyRingImpl::new(base, 7);

    let basis = ring.with_wrapped_indeterminates_dyn(|[X0, X1, X2, X3, X4, X5, X6]| [
        6 + 2 * X5 + 2 * X4 + X6 + 4 * X0 + 5 * X6 * X5 + X6 * X4 + 3 * X0 * X4 + 6 * X0 * X6 + 2 * X0 * X3 + X0 * X2 + 4 * X0 * X1 + 2 * X3 * X4 * X5 + 4 * X0 * X6 * X5 + 6 * X0 * X2 * X5 + 5 * X0 * X6 * X4 + 2 * X0 * X3 * X4 + 4 * X0 * X1 * X4 + X0 * X6.pow_ref(2) + 3 * X0 * X3 * X6 + 5 * X0 * X2 * X6 + 2 * X0 * X1 * X6 + X0 * X3.pow_ref(2) + 2 * X0 * X2 * X3 + 3 * X0 * X3 * X4 * X5 + 4 * X0 * X3 * X6 * X5 + 3 * X0 * X1 * X6 * X5 + 3 * X0 * X2 * X3 * X5 + 3 * X0 * X3 * X6 * X4 + 2 * X0 * X1 * X6 * X4 + 2 * X0 * X3.pow_ref(2) * X4 + 2 * X0 * X2 * X3 * X4 + 3 * X0 * X3.pow_ref(2) * X4 * X5 + 4 * X0 * X1 * X3 * X4 * X5 + X0 * X3.pow_ref(2) * X4.pow_ref(2),
        5 + 4 * X0 + 6 * X4 * X5 + 3 * X6 * X5 + 4 * X0 * X4 + 3 * X0 * X6 + 6 * X0 * X3 + 6 * X0 * X2 + 6 * X6 * X4 * X5 + 2 * X0 * X4 * X5 + 4 * X0 * X6 * X5 + 3 * X0 * X2 * X5 + 3 * X0 * X6 * X4 + 5 * X0 * X3 * X4 + 6 * X0 * X2 * X4 + 4 * X0 * X6.pow_ref(2) + 3 * X0 * X3 * X6 + 3 * X0 * X2 * X6 + 2 * X0 * X6 * X4 * X5 + 6 * X0 * X3 * X4 * X5 + 5 * X0 * X1 * X4 * X5 + 6 * X0 * X6.pow_ref(2) * X5 + 2 * X0 * X3 * X6 * X5 + 2 * X0 * X2 * X6 * X5 + 6 * X0 * X1 * X6 * X5 + 6 * X0 * X2 * X3 * X5 + 6 * X0 * X3 * X4.pow_ref(2) + 4 * X0 * X6.pow_ref(2) * X4 + 6 * X0 * X3 * X6 * X4 + 3 * X0 * X2 * X6 * X4 + 4 * X0 * X3 * X6 * X4 * X5 + 5 * X0 * X1 * X6 * X4 * X5 + 6 * X0 * X3.pow_ref(2) * X4 * X5 + 5 * X0 * X2 * X3 * X4 * X5 + 3 * X0 * X3 * X6 * X4.pow_ref(2) + 6 * X0 * X3.pow_ref(2) * X4.pow_ref(2) * X5.clone(),
        2 + 2 * X0 + 4 * X0 * X4 + 2 * X0 * X6 + 5 * X0 * X4 * X5 + 2 * X0 * X6 * X5 + 4 * X0 * X2 * X5 + 2 * X0 * X4.pow_ref(2) + 4 * X0 * X6 * X4 + 4 * X0 * X6.pow_ref(2) + 2 * X6 * X4 * X5.pow_ref(2) + 4 * X0 * X6 * X4 * X5 + X0 * X3 * X4 * X5 + X0 * X2 * X4 * X5 + 3 * X0 * X6.pow_ref(2) * X5 + 2 * X0 * X3 * X6 * X5 + 4 * X0 * X2 * X6 * X5 + 2 * X0 * X6 * X4.pow_ref(2) + X0 * X6.pow_ref(2) * X4 + 3 * X0 * X6 * X4 * X5.pow_ref(2) + 2 * X0 * X6.pow_ref(2) * X5.pow_ref(2) + 3 * X0 * X2 * X6 * X5.pow_ref(2) + X0 * X3 * X4.pow_ref(2) * X5 + X0 * X6.pow_ref(2) * X4 * X5 + X0 * X3 * X6 * X4 * X5 + 6 * X0 * X2 * X6 * X4 * X5 + 4 * X0 * X6.pow_ref(2) * X4.pow_ref(2) + 6 * X0 * X3 * X6 * X4 * X5.pow_ref(2) + 4 * X0 * X1 * X6 * X4 * X5.pow_ref(2) + 4 * X0 * X2 * X3 * X4 * X5.pow_ref(2) + 6 * X0 * X3 * X6 * X4.pow_ref(2) * X5 + 2 * X0 * X3.pow_ref(2) * X4.pow_ref(2) * X5.pow_ref(2),
        4 + 5 * X0 * X4 * X5 + 6 * X0 * X6 * X5 + 5 * X0 * X4.pow_ref(2) * X5 + 3 * X0 * X6 * X4 * X5 + 3 * X0 * X6.pow_ref(2) * X5 + 6 * X0 * X6 * X4 * X5.pow_ref(2) + 5 * X0 * X2 * X4 * X5.pow_ref(2) + X0 * X6.pow_ref(2) * X5.pow_ref(2) + 6 * X0 * X2 * X6 * X5.pow_ref(2) + 4 * X0 * X6 * X4.pow_ref(2) * X5 + 2 * X0 * X6.pow_ref(2) * X4 * X5 + 5 * X0 * X3 * X4.pow_ref(2) * X5.pow_ref(2) + 3 * X0 * X6.pow_ref(2) * X4 * X5.pow_ref(2) + 5 * X0 * X3 * X6 * X4 * X5.pow_ref(2) + 4 * X0 * X2 * X6 * X4 * X5.pow_ref(2) + 6 * X0 * X6.pow_ref(2) * X4.pow_ref(2) * X5 + 4 * X0 * X3 * X6 * X4.pow_ref(2) * X5.pow_ref(2),
        4 + 4 * X0 * X4.pow_ref(2) * X5.pow_ref(2) + X0 * X6 * X4 * X5.pow_ref(2) + X0 * X6.pow_ref(2) * X5.pow_ref(2) + 5 * X0 * X6 * X4.pow_ref(2) * X5.pow_ref(2) + 6 * X0 * X6.pow_ref(2) * X4 * X5.pow_ref(2) + 3 * X0 * X6.pow_ref(2) * X4 * X5.pow_ref(3) + 4 * X0 * X2 * X6 * X4 * X5.pow_ref(3) + 6 * X0 * X6.pow_ref(2) * X4.pow_ref(2) * X5.pow_ref(2) + 4 * X0 * X3 * X6 * X4.pow_ref(2) * X5.pow_ref(3),
        5 * X0 * X6 * X4.pow_ref(2) * X5.pow_ref(3) + 6 * X0 * X6.pow_ref(2) * X4 * X5.pow_ref(3) + 5 * X0 * X6.pow_ref(2) * X4.pow_ref(2) * X5.pow_ref(3),
        2 * X0 * X6.pow_ref(2) * X4.pow_ref(2) * X5.pow_ref(4)
    ]);

    let start = std::time::Instant::now();
    let gb = buchberger(&ring, basis, DegRevLex, default_sort_fn(&ring, DegRevLex), |_| false, TEST_COMPUTATION_CONTROLLER).unwrap_or_else(no_error);
    let end = std::time::Instant::now();

    println!("Computed GB in {} ms", (end - start).as_millis());

    assert_eq!(130, gb.len());
}

#[test]
#[ignore]
fn test_groebner_cyclic6() {
    let base = zn_static::Fp::<65537>::RING;
    let ring = MultivariatePolyRingImpl::new(base, 6);

    let cyclic6 = ring.with_wrapped_indeterminates_dyn(|[x, y, z, t, u, v]| {
        [x + y + z + t + u + v, x*y + y*z + z*t + t*u + x*v + u*v, x*y*z + y*z*t + z*t*u + x*y*v + x*u*v + t*u*v, x*y*z*t + y*z*t*u + x*y*z*v + x*y*u*v + x*t*u*v + z*t*u*v, x*y*z*t*u + x*y*z*t*v + x*y*z*u*v + x*y*t*u*v + x*z*t*u*v + y*z*t*u*v, x*y*z*t*u*v - 1]
    });

    let start = std::time::Instant::now();
    let gb = buchberger(&ring, cyclic6, DegRevLex, default_sort_fn(&ring, DegRevLex), |_| false, TEST_COMPUTATION_CONTROLLER).unwrap_or_else(no_error);
    let end = std::time::Instant::now();

    println!("Computed GB in {} ms", (end - start).as_millis());
    assert_eq!(45, gb.len());

}

#[test]
#[ignore]
fn test_groebner_cyclic7() {
    let base = zn_static::Fp::<65537>::RING;
    let ring = MultivariatePolyRingImpl::new(base, 7);

    let cyclic7 = ring.with_wrapped_indeterminates_dyn(|[x, y, z, t, u, v, w]| [
        x + y + z + t + u + v + w, x*y + y*z + z*t + t*u + u*v + x*w + v*w, x*y*z + y*z*t + z*t*u + t*u*v + x*y*w + x*v*w + u*v*w, x*y*z*t + y*z*t*u + z*t*u*v + x*y*z*w + x*y*v*w + x*u*v*w + t*u*v*w, 
        x*y*z*t*u + y*z*t*u*v + x*y*z*t*w + x*y*z*v*w + x*y*u*v*w + x*t*u*v*w + z*t*u*v*w, x*y*z*t*u*v + x*y*z*t*u*w + x*y*z*t*v*w + x*y*z*u*v*w + x*y*t*u*v*w + x*z*t*u*v*w + y*z*t*u*v*w, x*y*z*t*u*v*w - 1
    ]);

    let start = std::time::Instant::now();
    let gb = buchberger(&ring, cyclic7, DegRevLex, default_sort_fn(&ring, DegRevLex), |_| false, TEST_COMPUTATION_CONTROLLER).unwrap_or_else(no_error);
    let end = std::time::Instant::now();

    println!("Computed GB in {} ms", (end - start).as_millis());
    assert_eq!(209, gb.len());
}

#[test]
#[ignore]
fn test_groebner_cyclic8() {
    let base = zn_static::Fp::<65537>::RING;
    let ring = MultivariatePolyRingImpl::new(base, 8);

    let cyclic7 = ring.with_wrapped_indeterminates_dyn(|[x, y, z, s, t, u, v, w]| [
        x + y + z + s + t + u + v + w, x*y + y*z + z*s + s*t + t*u + u*v + x*w + v*w, x*y*z + y*z*s + z*s*t + s*t*u + t*u*v + x*y*w + x*v*w + u*v*w, 
        x*y*z*s + y*z*s*t + z*s*t*u + s*t*u*v + x*y*z*w + x*y*v*w + x*u*v*w + t*u*v*w, x*y*z*s*t + y*z*s*t*u + z*s*t*u*v + x*y*z*s*w + x*y*z*v*w + x*y*u*v*w + x*t*u*v*w + s*t*u*v*w, x*y*z*s*t*u + y*z*s*t*u*v + x*y*z*s*t*w + x*y*z*s*v*w + x*y*z*u*v*w + x*y*t*u*v*w + x*s*t*u*v*w + z*s*t*u*v*w, 
        x*y*z*s*t*u*v + x*y*z*s*t*u*w + x*y*z*s*t*v*w + x*y*z*s*u*v*w + x*y*z*t*u*v*w + x*y*s*t*u*v*w + x*z*s*t*u*v*w + y*z*s*t*u*v*w, x*y*z*s*t*u*v*w - 1
    ]);

    let start = std::time::Instant::now();
    let gb = buchberger(&ring, cyclic7, DegRevLex, default_sort_fn(&ring, DegRevLex), |_| false, TEST_COMPUTATION_CONTROLLER).unwrap_or_else(no_error);
    let end = std::time::Instant::now();

    println!("Computed GB in {} ms", (end - start).as_millis());
    assert_eq!(372, gb.len());
}