//! Multivariate polynomial structures and methods.
use ahash::{HashMap, HashMapExt};
use std::cell::{Cell, UnsafeCell};
use std::cmp::{Ordering, Reverse};
use std::collections::{BTreeMap, BinaryHeap};
use std::fmt::Display;
use std::marker::PhantomData;
use std::mem;
use std::ops::{Add, Div, Mul, Neg, Sub};
use std::sync::Arc;
use crate::domains::algebraic_number::AlgebraicExtension;
use crate::domains::float::FloatLike;
use crate::domains::integer::{Integer, IntegerRing};
use crate::domains::rational::{Fraction, FractionField, FractionNormalization, Q, RationalField};
use crate::domains::{
Derivable, EuclideanDomain, Field, InternalOrdering, Ring, RingOps, SelfRing, Set,
};
use crate::printer::{AtomPrinter, PrintOptions, PrintState};
use super::gcd::PolynomialGCD;
use super::univariate::UnivariatePolynomial;
use super::{Exponent, INLINED_EXPONENTS, LexOrder, MonomialOrder, PolyVariable, PositiveExponent};
use smallvec::{SmallVec, smallvec};
const MAX_DENSE_MUL_BUFFER_SIZE: usize = 1 << 24;
thread_local! { static DENSE_MUL_BUFFER: Cell<Vec<u32>> = const { Cell::new(Vec::new()) }; }
/// A ring for multivariate polynomials.
#[derive(Clone, PartialEq, Eq, Hash, Debug)]
pub struct PolynomialRing<R: Ring, E: Exponent = u16> {
pub(crate) ring: R,
_phantom_exp: PhantomData<E>,
}
impl<R: Ring + FractionNormalization, E: Exponent> FractionNormalization for PolynomialRing<R, E> {
fn get_normalization_factor(&self, a: &Self::Element) -> Self::Element {
a.constant(a.ring.get_normalization_factor(&a.lcoeff()))
}
}
impl<R: EuclideanDomain + FractionNormalization, E: Exponent> PolynomialRing<FractionField<R>, E> {
pub fn to_rational_polynomial(
&self,
e: &<Self as Set>::Element,
) -> Fraction<PolynomialRing<R, E>> {
let ring = self.ring.ring();
let mut lcm = ring.one();
for x in &e.coefficients {
let g = ring.gcd(&lcm, x.denominator_ref());
lcm = ring.mul(&lcm, &ring.quot_rem(x.denominator_ref(), &g).0);
}
let e2 = e.map_coeff(
|c| {
ring.mul(
c.numerator_ref(),
&ring.quot_rem(&lcm, c.denominator_ref()).0,
)
},
ring.clone(),
);
Fraction::from_unchecked(e2.constant(lcm), e2)
}
}
impl<R: Ring, E: Exponent> PolynomialRing<R, E> {
pub fn new(coeff_ring: R) -> PolynomialRing<R, E> {
PolynomialRing {
ring: coeff_ring,
_phantom_exp: PhantomData,
}
}
pub fn from_poly(poly: &MultivariatePolynomial<R, E>) -> PolynomialRing<R, E> {
PolynomialRing {
ring: poly.ring.clone(),
_phantom_exp: PhantomData,
}
}
/// Get the coefficient ring.
pub fn coefficient_ring(&self) -> &R {
&self.ring
}
}
impl<R: Ring, E: Exponent> std::fmt::Display for PolynomialRing<R, E> {
fn fmt(&self, _: &mut std::fmt::Formatter<'_>) -> std::fmt::Result {
Ok(())
}
}
impl<R: Ring, E: Exponent> Set for PolynomialRing<R, E> {
type Element = MultivariatePolynomial<R, E>;
fn size(&self) -> Option<Integer> {
None
}
}
impl<R: Ring, E: Exponent> RingOps<MultivariatePolynomial<R, E>> for PolynomialRing<R, E> {
#[inline]
fn add(&self, a: Self::Element, b: Self::Element) -> Self::Element {
a + b
}
#[inline]
fn sub(&self, a: Self::Element, b: Self::Element) -> Self::Element {
a - b
}
#[inline]
fn mul(&self, a: Self::Element, b: Self::Element) -> Self::Element {
a * &b
}
#[inline]
fn add_assign(&self, a: &mut Self::Element, b: Self::Element) {
*a = std::mem::replace(a, b.zero()) + b;
}
#[inline]
fn sub_assign(&self, a: &mut Self::Element, b: Self::Element) {
*a = std::mem::replace(a, b.zero()) - b;
}
#[inline]
fn mul_assign(&self, a: &mut Self::Element, b: Self::Element) {
*a = &*a * &b;
}
#[inline]
fn add_mul_assign(&self, a: &mut Self::Element, b: Self::Element, c: Self::Element) {
*a = std::mem::replace(a, b.zero()) + b * &c
}
#[inline]
fn sub_mul_assign(&self, a: &mut Self::Element, b: Self::Element, c: Self::Element) {
*a = std::mem::replace(a, b.zero()) - b * &c
}
#[inline]
fn neg(&self, a: Self::Element) -> Self::Element {
a.neg()
}
}
impl<R: Ring, E: Exponent> RingOps<&MultivariatePolynomial<R, E>> for PolynomialRing<R, E> {
#[inline]
fn add(&self, a: &Self::Element, b: &Self::Element) -> Self::Element {
a + b
}
#[inline]
fn sub(&self, a: &Self::Element, b: &Self::Element) -> Self::Element {
a - b
}
#[inline]
fn mul(&self, a: &Self::Element, b: &Self::Element) -> Self::Element {
a * b
}
#[inline]
fn add_assign(&self, a: &mut Self::Element, b: &Self::Element) {
*a = &*a + b;
}
#[inline]
fn sub_assign(&self, a: &mut Self::Element, b: &Self::Element) {
*a = &*a - b;
}
#[inline]
fn mul_assign(&self, a: &mut Self::Element, b: &Self::Element) {
*a = std::mem::replace(a, b.zero()) * b;
}
#[inline]
fn add_mul_assign(&self, a: &mut Self::Element, b: &Self::Element, c: &Self::Element) {
*a = std::mem::replace(a, b.zero()) + b * c
}
#[inline]
fn sub_mul_assign(&self, a: &mut Self::Element, b: &Self::Element, c: &Self::Element) {
*a = std::mem::replace(a, b.zero()) - b * c
}
#[inline]
fn neg(&self, a: &Self::Element) -> Self::Element {
a.clone().neg()
}
}
impl<R: Ring, E: Exponent> Ring for PolynomialRing<R, E> {
#[inline]
fn zero(&self) -> Self::Element {
MultivariatePolynomial::new(&self.ring, None, Arc::new(vec![]))
}
#[inline]
fn one(&self) -> Self::Element {
self.zero().one()
}
#[inline]
fn nth(&self, n: Integer) -> Self::Element {
self.zero().constant(self.ring.nth(n))
}
#[inline]
fn pow(&self, b: &Self::Element, e: u64) -> Self::Element {
b.pow(e as usize)
}
#[inline]
fn is_zero(&self, a: &Self::Element) -> bool {
a.is_zero()
}
#[inline]
fn is_one(&self, a: &Self::Element) -> bool {
a.is_one()
}
fn one_is_gcd_unit() -> bool {
false
}
fn characteristic(&self) -> Integer {
self.ring.characteristic()
}
fn try_inv(&self, a: &Self::Element) -> Option<Self::Element> {
if a.is_constant() {
let inv = self.ring.try_inv(&a.get_constant())?;
Some(a.constant(inv))
} else {
None
}
}
fn try_div(&self, a: &Self::Element, b: &Self::Element) -> Option<Self::Element> {
a.try_div(b)
}
fn sample(&self, _rng: &mut impl rand::RngCore, _range: (i64, i64)) -> Self::Element {
todo!("Sampling a polynomial is not possible yet")
}
fn format<W: std::fmt::Write>(
&self,
element: &Self::Element,
opts: &PrintOptions,
state: PrintState,
f: &mut W,
) -> Result<bool, std::fmt::Error> {
element.format(opts, state, f)
}
fn has_independent_elements(&self) -> bool {
// the coefficient ring is stored in the polynomial
true
}
}
impl<R: EuclideanDomain + PolynomialGCD<E>, E: PositiveExponent> EuclideanDomain
for PolynomialRing<R, E>
{
fn rem(&self, a: &Self::Element, b: &Self::Element) -> Self::Element {
a.rem(b)
}
fn quot_rem(&self, a: &Self::Element, b: &Self::Element) -> (Self::Element, Self::Element) {
a.quot_rem(b, false)
}
fn gcd(&self, a: &Self::Element, b: &Self::Element) -> Self::Element {
a.gcd(b)
}
}
/// Multivariate polynomial with a sparse degree and dense variable representation.
/// Negative exponents are supported, if they are allowed by the exponent type.
#[derive(Clone)]
#[cfg_attr(feature = "bincode", derive(bincode_trait_derive::Encode))]
pub struct MultivariatePolynomial<F: Ring, E: Exponent = u16, O: MonomialOrder = LexOrder> {
// Data format: the i-th monomial is stored as coefficients[i] and
// exponents[i * nvars .. (i + 1) * nvars]. Terms are always expanded and sorted by the exponents via
// cmp_exponents().
pub coefficients: Vec<F::Element>,
pub exponents: Vec<E>,
/// The coefficient ring.
pub ring: F,
pub variables: Arc<Vec<PolyVariable>>,
pub(crate) _phantom: PhantomData<O>,
}
#[cfg(feature = "bincode")]
impl<
C: crate::state::HasStateMap,
F: Ring + bincode::Decode<C>,
E: Exponent + bincode::Decode<C>,
O: MonomialOrder + bincode::Decode<C>,
> bincode::Decode<C> for MultivariatePolynomial<F, E, O>
where
F::Element: for<'a> bincode::Decode<&'a F>,
{
fn decode<D: bincode::de::Decoder<Context = C>>(
decoder: &mut D,
) -> Result<Self, bincode::error::DecodeError> {
let ring = F::decode(decoder)?;
let coefficients = Vec::<F::Element>::decode(&mut decoder.with_context(&ring))?;
let exponents = Vec::<E>::decode(decoder)?;
let variables = Arc::<Vec<PolyVariable>>::decode(decoder)?;
Ok(MultivariatePolynomial {
coefficients,
exponents,
ring,
variables,
_phantom: PhantomData,
})
}
}
impl<F: Ring, E: Exponent, O: MonomialOrder> MultivariatePolynomial<F, E, O> {
/// Constructs a zero polynomial. Instead of using this constructor,
/// prefer to create new polynomials from existing ones, so that the
/// variable map and field are inherited.
#[inline]
pub fn new(ring: &F, cap: Option<usize>, variables: Arc<Vec<PolyVariable>>) -> Self {
Self {
coefficients: Vec::with_capacity(cap.unwrap_or(0)),
exponents: Vec::with_capacity(cap.unwrap_or(0) * variables.len()),
ring: ring.clone(),
variables,
_phantom: PhantomData,
}
}
/// Constructs a zero polynomial. Instead of using this constructor,
/// prefer to create new polynomials from existing ones, so that the
/// variable map is inherited.
#[inline]
pub fn new_zero(ring: &F) -> Self {
Self {
coefficients: vec![],
exponents: vec![],
ring: ring.clone(),
variables: Arc::new(vec![]),
_phantom: PhantomData,
}
}
/// Constructs a polynomial that is one. Instead of using this constructor,
/// prefer to create new polynomials from existing ones, so that the
/// variable map is inherited.
#[inline]
pub fn new_one(ring: &F) -> Self {
Self {
coefficients: vec![ring.one()],
exponents: vec![],
ring: ring.clone(),
variables: Arc::new(vec![]),
_phantom: PhantomData,
}
}
/// Constructs a zero polynomial, inheriting the field and variable map from `self`.
#[inline]
pub fn zero(&self) -> Self {
Self {
coefficients: vec![],
exponents: vec![],
ring: self.ring.clone(),
variables: self.variables.clone(),
_phantom: PhantomData,
}
}
/// Constructs a zero polynomial with the given number of variables and capacity,
/// inheriting the field and variable map from `self`.
#[inline]
pub fn zero_with_capacity(&self, cap: usize) -> Self {
Self {
coefficients: Vec::with_capacity(cap),
exponents: Vec::with_capacity(cap * self.nvars()),
ring: self.ring.clone(),
variables: self.variables.clone(),
_phantom: PhantomData,
}
}
/// Constructs a constant polynomial,
/// inheriting the field and variable map from `self`.
#[inline]
pub fn constant(&self, coeff: F::Element) -> Self {
if self.ring.is_zero(&coeff) {
return self.zero();
}
Self {
coefficients: vec![coeff],
exponents: vec![E::zero(); self.nvars()],
ring: self.ring.clone(),
variables: self.variables.clone(),
_phantom: PhantomData,
}
}
/// Constructs a polynomial that is one, inheriting the field and variable map from `self`.
#[inline]
pub fn one(&self) -> Self {
Self {
coefficients: vec![self.ring.one()],
exponents: vec![E::zero(); self.nvars()],
ring: self.ring.clone(),
variables: self.variables.clone(),
_phantom: PhantomData,
}
}
/// Constructs a polynomial with a single term.
#[inline]
pub fn monomial(&self, coeff: F::Element, exponents: Vec<E>) -> Self {
debug_assert!(self.nvars() == exponents.len());
if self.ring.is_zero(&coeff) {
return self.zero();
}
Self {
coefficients: vec![coeff],
exponents,
ring: self.ring.clone(),
variables: self.variables.clone(),
_phantom: PhantomData,
}
}
/// Constructs a polynomial with a single term that is a variable.
#[inline]
pub fn variable(&self, var: &PolyVariable) -> Result<Self, String> {
if let Some(pos) = self.variables.iter().position(|v| v == var) {
let mut exp = vec![E::zero(); self.nvars()];
exp[pos] = E::one();
Ok(self.monomial(self.ring.one(), exp))
} else {
Err(format!("Variable {} not found", var))
}
}
/// Get the ith monomial
pub fn to_monomial_view(&self, i: usize) -> MonomialView<'_, F, E> {
assert!(i < self.nterms());
MonomialView {
coefficient: &self.coefficients[i],
exponents: self.exponents(i),
}
}
#[inline]
pub fn reserve(&mut self, cap: usize) -> &mut Self {
self.coefficients.reserve(cap);
self.exponents.reserve(cap * self.nvars());
self
}
#[inline]
pub fn is_zero(&self) -> bool {
self.nterms() == 0
}
#[inline]
pub fn is_one(&self) -> bool {
self.nterms() == 1
&& self.ring.is_one(&self.coefficients[0])
&& self.exponents.iter().all(|x| x.is_zero())
}
/// Returns the number of terms in the polynomial.
#[inline]
pub fn nterms(&self) -> usize {
self.coefficients.len()
}
/// Returns the number of variables in the polynomial.
#[inline]
pub fn nvars(&self) -> usize {
self.variables.len()
}
/// Returns true if the polynomial is constant.
#[inline]
pub fn is_constant(&self) -> bool {
if self.is_zero() {
return true;
}
if self.nterms() >= 2 {
return false;
}
debug_assert!(!self.ring.is_zero(self.coefficients.first().unwrap()));
self.exponents.iter().all(|e| e.is_zero())
}
/// Get the constant term of the polynomial.
#[inline]
pub fn get_constant(&self) -> F::Element {
if self.is_zero() || !self.exponents(0).iter().all(|e| e.is_zero()) {
return self.ring.zero();
}
self.coefficients[0].clone()
}
/// Returns the `index`th monomial, starting from the back.
#[inline]
pub fn coefficient_back(&self, index: usize) -> &F::Element {
&self.coefficients[self.nterms() - index - 1]
}
/// Returns the slice for the exponents of the specified monomial.
#[inline]
pub fn exponents(&self, index: usize) -> &[E] {
//&self.exponents[index * self.nvars()..(index + 1) * self.nvars()]
unsafe {
self.exponents
.get_unchecked(index * self.nvars()..(index + 1) * self.nvars())
}
}
/// Returns the slice for the exponents of the specified monomial
/// starting from the back.
#[inline]
pub fn exponents_back(&self, index: usize) -> &[E] {
let index = self.nterms() - index - 1;
&self.exponents[index * self.nvars()..(index + 1) * self.nvars()]
}
#[inline(always)]
pub fn last_exponents(&self) -> &[E] {
//assert!(self.nterms() > 0);
&self.exponents[(self.nterms() - 1) * self.nvars()..self.nterms() * self.nvars()]
}
/// Returns the mutable slice for the exponents of the specified monomial.
#[inline]
pub fn exponents_mut(&mut self, index: usize) -> &mut [E] {
let nvars = self.nvars();
&mut self.exponents[index * nvars..(index + 1) * nvars]
}
/// Returns an iterator over the exponents of every monomial.
#[inline]
pub fn exponents_iter(&self) -> std::slice::Chunks<'_, E> {
self.exponents.chunks(self.nvars())
}
/// Returns an iterator over the mutable exponents of every monomial.
#[inline]
pub fn exponents_iter_mut(&mut self) -> std::slice::ChunksMut<'_, E> {
let nvars = self.nvars();
self.exponents.chunks_mut(nvars)
}
/// Reset the polynomial to 0.
#[inline]
pub fn clear(&mut self) {
self.coefficients.clear();
self.exponents.clear();
}
/// Get a copy of the variable list.
pub fn get_vars(&self) -> Arc<Vec<PolyVariable>> {
self.variables.clone()
}
/// Get a reference to the variables list.
pub fn get_vars_ref(&self) -> &[PolyVariable] {
self.variables.as_ref()
}
/// Rename a variable.
pub fn rename_variable(&mut self, old: &PolyVariable, new: &PolyVariable) {
if let Some(pos) = self.variables.iter().position(|v| v == old) {
let mut new_vars = self.variables.as_ref().clone();
new_vars[pos] = new.clone();
self.variables = Arc::new(new_vars);
}
}
/// Unify the variable maps of two polynomials, i.e.
/// rewrite a polynomial in `x` and one in `y` to a
/// two polynomial in `x` and `y`.
///
/// The variable map will be inherited from
/// `self` and will be extended by variables occurring
/// in `other`.
#[inline(always)]
pub fn unify_variables(&mut self, other: &mut Self) {
if self.variables == other.variables {
// reduce the number of equal variable copies
if self.variables.as_ptr() != other.variables.as_ptr() {
if self.variables.as_ptr() < other.variables.as_ptr() {
other.variables = self.variables.clone();
} else {
self.variables = other.variables.clone();
}
}
return;
}
self.unify_variables_impl(other)
}
fn unify_variables_impl(&mut self, other: &mut Self) {
let mut new_var_map = self.variables.as_ref().clone();
let mut new_var_pos_other = vec![0; other.nvars()];
for (pos, v) in new_var_pos_other.iter_mut().zip(other.variables.as_ref()) {
if let Some(p) = new_var_map.iter().position(|x| x == v) {
*pos = p;
} else {
*pos = new_var_map.len();
new_var_map.push(v.clone());
}
}
let mut newexp = vec![E::zero(); new_var_map.len() * self.nterms()];
for t in 0..self.nterms() {
newexp[t * new_var_map.len()..t * new_var_map.len() + self.nvars()]
.copy_from_slice(self.exponents(t));
}
self.variables = Arc::new(new_var_map);
self.exponents = newexp;
// check if term ordering remains unchanged
if new_var_pos_other.windows(2).all(|w| w[0] <= w[1]) {
let mut newexp = vec![E::zero(); self.nvars() * other.nterms()];
if other.nvars() > 0 {
for (d, t) in newexp
.chunks_mut(self.nvars())
.zip(other.exponents.chunks(other.nvars()))
{
for (var, e) in t.iter().enumerate() {
d[new_var_pos_other[var]] = *e;
}
}
}
other.variables = self.variables.clone();
other.exponents = newexp;
return;
}
// reconstruct 'other' with correct monomial ordering
let mut newother = Self::new(&other.ring, other.nterms().into(), self.variables.clone());
let mut newexp = vec![E::zero(); self.nvars()];
for t in other.into_iter() {
for c in &mut newexp {
*c = E::zero();
}
for (var, e) in t.exponents.iter().enumerate() {
newexp[new_var_pos_other[var]] = *e;
}
newother.append_monomial(t.coefficient.clone(), &newexp);
}
*other = newother;
}
/// Unify the variable maps of all polynomials in the slice.
pub fn unify_variables_list(polys: &mut [Self]) {
if polys.len() < 2 {
return;
}
let (first, rest) = polys.split_first_mut().unwrap();
for _ in 0..2 {
for p in &mut *rest {
first.unify_variables(p);
}
}
}
/// Reverse the coefficients: `1+2x+3x^3` becomes `3+2x^2+x^3`.
pub fn reverse(&mut self) {
self.coefficients.reverse();
let nterms = self.nterms();
let nvars = self.nvars();
let degs = (0..nvars).map(|i| self.degree(i)).collect::<Vec<_>>();
for e in self.exponents.chunks_mut(nvars) {
for (ee, d) in e.iter_mut().zip(°s) {
*ee = *d - *ee;
}
}
let midu = if nterms.is_multiple_of(2) {
self.nvars() * (nterms / 2)
} else {
self.nvars() * (nterms / 2 + 1)
};
let (l, r) = self.exponents.split_at_mut(midu);
let rend = r.len();
for i in 0..nterms / 2 {
l[i * nvars..(i + 1) * nvars]
.swap_with_slice(&mut r[rend - (i + 1) * nvars..rend - i * nvars]);
}
}
/// Reverse the monomial ordering in-place.
fn reverse_monomials(&mut self) {
let nterms = self.nterms();
let nvars = self.nvars();
if nterms < 2 {
return;
}
self.coefficients.reverse();
let midu = if nterms.is_multiple_of(2) {
self.nvars() * (nterms / 2)
} else {
self.nvars() * (nterms / 2 + 1)
};
let (l, r) = self.exponents.split_at_mut(midu);
let rend = r.len();
for i in 0..nterms / 2 {
l[i * nvars..(i + 1) * nvars]
.swap_with_slice(&mut r[rend - (i + 1) * nvars..rend - i * nvars]);
}
}
/// Add a variable to the polynomial if it is not already present.
pub fn add_variable(&mut self, var: &PolyVariable) {
if self.variables.iter().any(|v| v == var) {
return;
}
let l = self.variables.len();
let mut new_exp = vec![E::zero(); (l + 1) * self.nterms()];
if l > 0 {
for (en, e) in new_exp.chunks_mut(l + 1).zip(self.exponents.chunks(l)) {
en[..l].copy_from_slice(e);
}
}
let mut new_vars = self.variables.as_ref().clone();
new_vars.push(var.clone());
self.variables = Arc::new(new_vars);
self.exponents = new_exp;
}
/// Add several variables to the polynomial if they are not already present.
pub fn add_variables(&mut self, vars: &[PolyVariable]) {
// collect only genuinely new variables
let new_vars: Vec<_> = vars
.iter()
.filter(|var| !self.variables.iter().any(|v| v == *var))
.cloned()
.collect();
if new_vars.is_empty() {
return;
}
let l_old = self.variables.len();
let n_new = new_vars.len();
let l_new = l_old + n_new;
let mut new_exp = vec![E::zero(); l_new * self.nterms()];
if l_old > 0 {
for (en, e) in new_exp.chunks_mut(l_new).zip(self.exponents.chunks(l_old)) {
en[..l_old].copy_from_slice(e);
}
}
let mut variables = self.variables.as_ref().clone();
variables.extend(new_vars);
self.variables = Arc::new(variables);
self.exponents = new_exp;
}
/// Check if the polynomial is sorted and has only non-zero coefficients
pub fn check_consistency(&self) {
assert_eq!(self.coefficients.len(), self.nterms());
assert_eq!(self.exponents.len(), self.nterms() * self.nvars());
for c in &self.coefficients {
if self.ring.is_zero(c) {
panic!("Inconsistent polynomial (0 coefficient): {self}");
}
}
for t in 1..self.nterms() {
match O::cmp(self.exponents(t), self.exponents(t - 1)) {
Ordering::Equal => panic!("Inconsistent polynomial (equal monomials): {self}"),
Ordering::Less => {
panic!("Inconsistent polynomial (wrong monomial ordering): {self}")
}
Ordering::Greater => {}
}
}
}
/// Append a monomial to the back. It merges with the last monomial if the
/// exponents are equal.
#[inline]
pub fn append_monomial_back(&mut self, coefficient: F::Element, exponents: &[E]) {
if self.ring.is_zero(&coefficient) {
return;
}
let nterms = self.nterms();
if nterms > 0 && exponents == self.last_exponents() {
self.ring
.add_assign(&mut self.coefficients[nterms - 1], &coefficient);
if self.ring.is_zero(&self.coefficients[nterms - 1]) {
self.coefficients.pop();
self.exponents.truncate((nterms - 1) * self.nvars());
}
} else {
self.coefficients.push(coefficient);
self.exponents.extend_from_slice(exponents);
}
}
/// Appends a monomial to the polynomial.
pub fn append_monomial(&mut self, coefficient: F::Element, exponents: &[E]) {
if self.ring.is_zero(&coefficient) {
return;
}
if self.nvars() != exponents.len() {
panic!(
"nvars mismatched: got {}, expected {}",
exponents.len(),
self.nvars()
);
}
// should we append to the back?
if self.nterms() == 0 || O::cmp(self.last_exponents(), exponents).is_lt() {
self.coefficients.push(coefficient);
self.exponents.extend_from_slice(exponents);
return;
}
if O::cmp(self.exponents(0), exponents).is_gt() {
self.coefficients.insert(0, coefficient);
self.exponents.splice(0..0, exponents.iter().cloned());
return;
}
// Binary search to find the insert-point.
let mut l = 0;
let mut r = self.nterms();
while l <= r {
let m = (l + r) / 2;
let c = O::cmp(exponents, self.exponents(m)); // note the reversal
match c {
Ordering::Equal => {
// Add the two coefficients.
self.ring
.add_assign(&mut self.coefficients[m], &coefficient);
if self.ring.is_zero(&self.coefficients[m]) {
// The coefficient becomes zero. Remove this monomial.
self.coefficients.remove(m);
let i = m * self.nvars();
self.exponents.splice(i..i + self.nvars(), Vec::new());
}
return;
}
Ordering::Greater => {
l = m + 1;
if l == self.nterms() {
self.coefficients.push(coefficient);
self.exponents.extend_from_slice(exponents);
return;
}
}
Ordering::Less => {
if m == 0 {
self.coefficients.insert(0, coefficient);
self.exponents.splice(0..0, exponents.iter().cloned());
return;
}
r = m - 1;
}
}
}
self.coefficients.insert(l, coefficient);
let i = l * self.nvars();
self.exponents.splice(i..i, exponents.iter().cloned());
}
}
impl<F: Ring, E: Exponent, O: MonomialOrder> SelfRing for MultivariatePolynomial<F, E, O> {
#[inline]
fn is_zero(&self) -> bool {
self.is_zero()
}
#[inline]
fn is_one(&self) -> bool {
self.is_one()
}
fn format<W: std::fmt::Write>(
&self,
opts: &PrintOptions,
mut state: PrintState,
f: &mut W,
) -> Result<bool, std::fmt::Error> {
let print_ring = opts.print_ring && !self.ring.has_independent_elements();
if self.is_constant() {
if self.is_zero() {
if state.in_sum {
f.write_str("+")?;
}
f.write_char('0')?;
return Ok(false);
} else if !print_ring || state.level > 0 {
return self.ring.format(&self.coefficients[0], opts, state, f);
}
}
let add_paren = (self.nterms() > 1 || print_ring) && state.in_product
|| ((state.in_exp || state.in_exp_base)
&& (self.nterms() > 1
|| print_ring
|| self.exponents(0).iter().filter(|e| **e > E::zero()).count() > 1
|| !self.ring.is_one(&self.coefficients[0])));
if add_paren {
if state.in_sum {
f.write_str("+")?;
state.in_sum = false;
}
state.in_product = false;
state.in_exp = false;
state.in_exp_base = false;
f.write_str("(")?;
}
let in_product = state.in_product;
let var_map: Vec<String> = self
.variables
.as_ref()
.iter()
.map(|v| {
v.format_string(
opts,
PrintState {
in_exp: true,
..state
},
)
})
.collect();
state.level += 1;
for monomial in self {
let has_var = monomial.exponents.iter().any(|e| !e.is_zero());
state.in_product = in_product || has_var;
state.suppress_one = has_var; // any products before should not be considered
state.in_exp |= print_ring; // make sure to add parentheses
let mut suppressed_one = self.ring.format(monomial.coefficient, opts, state, f)?;
for (var_id, e) in var_map.iter().zip(monomial.exponents) {
if e.is_zero() {
continue;
}
if suppressed_one {
suppressed_one = false;
} else if !opts.mode.is_latex() {
f.write_char(opts.multiplication_operator)?;
}
f.write_str(var_id)?;
if e.to_i32() != 1 {
if opts.mode.is_latex() {
write!(f, "^{{{e}}}")?;
} else if opts.double_star_for_exponentiation {
write!(f, "**{e}")?;
} else if opts.mode.is_symbolica() && opts.num_exp_as_superscript {
state.superscript = true;
AtomPrinter::format_digits(e.to_string(), opts, &state, f)?;
state.superscript = false;
} else {
write!(f, "^{e}")?;
}
}
}
state.in_sum = true;
}
state.level -= 1;
if self.is_zero() {
f.write_char('0')?;
}
if print_ring && state.level == 0 {
self.ring.format_ring(opts, state, f)?;
}
if add_paren {
f.write_str(")")?;
}
Ok(false)
}
}
impl<F: Ring + std::fmt::Debug, E: Exponent + std::fmt::Debug, O: MonomialOrder> std::fmt::Debug
for MultivariatePolynomial<F, E, O>
{
fn fmt(&self, f: &mut std::fmt::Formatter) -> std::fmt::Result {
if self.is_zero() {
return write!(f, "[]");
}
let mut first = true;
write!(f, "[ ")?;
for monomial in self {
if first {
first = false;
} else {
write!(f, ", ")?;
}
write!(
f,
"{{ {:?}, {:?} }}",
monomial.coefficient, monomial.exponents
)?;
}
write!(f, " ]")
}
}
impl<F: Ring + Display, E: Exponent, O: MonomialOrder> Display for MultivariatePolynomial<F, E, O> {
fn fmt(&self, f: &mut std::fmt::Formatter) -> std::fmt::Result {
self.format(&PrintOptions::from_fmt(f), PrintState::from_fmt(f), f)
.map(|_| ())
}
}
impl<F: Ring, E: Exponent, O: MonomialOrder> PartialEq for MultivariatePolynomial<F, E, O> {
#[inline]
fn eq(&self, other: &Self) -> bool {
if self.variables != other.variables {
if self.is_constant() != other.is_constant() {
return false;
}
if self.is_zero() != other.is_zero() {
return false;
}
if self.is_zero() {
return true;
}
if self.is_constant() {
return self.coefficients[0] == other.coefficients[0];
}
return false;
}
if self.nterms() != other.nterms() {
return false;
}
self.exponents.eq(&other.exponents) && self.coefficients.eq(&other.coefficients)
}
}
impl<F: Ring, E: Exponent, O: MonomialOrder> std::hash::Hash for MultivariatePolynomial<F, E, O> {
fn hash<H: std::hash::Hasher>(&self, state: &mut H) {
self.coefficients.hash(state);
self.exponents.hash(state);
if !self.is_constant() {
self.variables.hash(state);
}
}
}
impl<F: Ring, E: Exponent, O: MonomialOrder> Eq for MultivariatePolynomial<F, E, O> {}
impl<R: Ring, E: Exponent, O: MonomialOrder> InternalOrdering for MultivariatePolynomial<R, E, O> {
/// An ordering of polynomials that has no intuitive meaning.
fn internal_cmp(&self, other: &Self) -> Ordering {
// TODO: what about different variables?
Ord::cmp(&self.exponents, &other.exponents)
.then_with(|| self.coefficients.internal_cmp(&other.coefficients))
}
}
impl<F: Ring, E: Exponent, O: MonomialOrder> Add for MultivariatePolynomial<F, E, O> {
type Output = Self;
fn add(mut self, mut other: Self) -> Self::Output {
assert_eq!(self.ring, other.ring);
self.unify_variables(&mut other);
if self.is_zero() {
return other;
}
if other.is_zero() {
return self;
}
// Merge the two polynomials, which are assumed to be already sorted.
let mut new_coefficients = vec![self.ring.zero(); self.nterms() + other.nterms()];
let mut new_exponents: Vec<E> =
vec![E::zero(); self.nvars() * (self.nterms() + other.nterms())];
let mut new_nterms = 0;
let mut i = 0;
let mut j = 0;
macro_rules! insert_monomial {
($source:expr, $index:expr) => {
mem::swap(
&mut new_coefficients[new_nterms],
&mut $source.coefficients[$index],
);
new_exponents[new_nterms * $source.nvars()..(new_nterms + 1) * $source.nvars()]
.clone_from_slice($source.exponents($index));
new_nterms += 1;
};
}
while i < self.nterms() && j < other.nterms() {
let c = O::cmp(self.exponents(i), other.exponents(j));
match c {
Ordering::Less => {
insert_monomial!(self, i);
i += 1;
}
Ordering::Greater => {
insert_monomial!(other, j);
j += 1;
}
Ordering::Equal => {
self.ring
.add_assign(&mut self.coefficients[i], &other.coefficients[j]);
if !self.ring.is_zero(&self.coefficients[i]) {
insert_monomial!(self, i);
}
i += 1;
j += 1;
}
}
}
while i < self.nterms() {
insert_monomial!(self, i);
i += 1;
}
while j < other.nterms() {
insert_monomial!(other, j);
j += 1;
}
new_coefficients.truncate(new_nterms);
new_exponents.truncate(self.nvars() * new_nterms);
Self {
coefficients: new_coefficients,
exponents: new_exponents,
ring: self.ring,
variables: self.variables,
_phantom: PhantomData,
}
}
}
impl<F: Ring, E: Exponent, O: MonomialOrder> Add<&MultivariatePolynomial<F, E, O>>
for &MultivariatePolynomial<F, E, O>
{
type Output = MultivariatePolynomial<F, E, O>;
fn add(self, other: &MultivariatePolynomial<F, E, O>) -> Self::Output {
assert_eq!(self.ring, other.ring);
if self.is_zero() {
return other.clone();
}
if other.is_zero() {
return self.clone();
}
if self.variables != other.variables {
let mut c1 = self.clone();
let mut c2 = other.clone();
c1.unify_variables(&mut c2);
return c1 + c2;
}
// Merge the two polynomials, which are assumed to be already sorted.
let mut new_coefficients = vec![self.ring.zero(); self.nterms() + other.nterms()];
let mut new_exponents: Vec<E> =
vec![E::zero(); self.nvars() * (self.nterms() + other.nterms())];
let mut new_nterms = 0;
let mut i = 0;
let mut j = 0;
macro_rules! insert_monomial {
($source:expr, $index:expr) => {
new_coefficients[new_nterms] = $source.coefficients[$index].clone();
new_exponents[new_nterms * $source.nvars()..(new_nterms + 1) * $source.nvars()]
.clone_from_slice($source.exponents($index));
new_nterms += 1;
};
}
while i < self.nterms() && j < other.nterms() {
let c = O::cmp(self.exponents(i), other.exponents(j));
match c {
Ordering::Less => {
insert_monomial!(self, i);
i += 1;
}
Ordering::Greater => {
insert_monomial!(other, j);
j += 1;
}
Ordering::Equal => {
let coeff = self.ring.add(&self.coefficients[i], &other.coefficients[j]);
if !self.ring.is_zero(&coeff) {
new_coefficients[new_nterms] = coeff;
new_exponents[new_nterms * self.nvars()..(new_nterms + 1) * self.nvars()]
.clone_from_slice(self.exponents(i));
new_nterms += 1;
}
i += 1;
j += 1;
}
}
}
while i < self.nterms() {
insert_monomial!(self, i);
i += 1;
}
while j < other.nterms() {
insert_monomial!(other, j);
j += 1;
}
new_coefficients.truncate(new_nterms);
new_exponents.truncate(self.nvars() * new_nterms);
MultivariatePolynomial {
coefficients: new_coefficients,
exponents: new_exponents,
ring: self.ring.clone(),
variables: self.variables.clone(),
_phantom: PhantomData,
}
}
}
impl<F: Ring, E: Exponent, O: MonomialOrder> Sub for MultivariatePolynomial<F, E, O> {
type Output = Self;
fn sub(self, other: Self) -> Self::Output {
self.add(other.neg())
}
}
impl<'a, F: Ring, E: Exponent, O: MonomialOrder> Sub<&'a MultivariatePolynomial<F, E, O>>
for &MultivariatePolynomial<F, E, O>
{
type Output = MultivariatePolynomial<F, E, O>;
fn sub(self, other: &'a MultivariatePolynomial<F, E, O>) -> Self::Output {
self + &other.clone().neg() // TODO: improve
}
}
impl<F: Ring, E: Exponent, O: MonomialOrder> Neg for MultivariatePolynomial<F, E, O> {
type Output = Self;
fn neg(mut self) -> Self::Output {
// Negate coefficients of all terms.
for c in &mut self.coefficients {
*c = self.ring.neg(&*c);
}
self
}
}
impl<'a, F: Ring, E: Exponent> Mul<&'a MultivariatePolynomial<F, E, LexOrder>>
for &MultivariatePolynomial<F, E, LexOrder>
{
type Output = MultivariatePolynomial<F, E, LexOrder>;
#[inline]
fn mul(self, rhs: &'a MultivariatePolynomial<F, E, LexOrder>) -> Self::Output {
assert_eq!(self.ring, rhs.ring);
if self.nterms() == 0 || rhs.nterms() == 0 {
return self.zero();
}
if self.is_constant() {
return rhs.clone().mul_coeff(self.coefficients[0].clone());
}
if rhs.is_constant() {
return self.clone().mul_coeff(rhs.coefficients[0].clone());
}
if self.variables != rhs.variables {
let mut c1 = self.clone();
let mut c2 = rhs.clone();
c1.unify_variables(&mut c2);
return c1.mul(&c2);
}
if self.nterms() == 1 {
return rhs
.clone()
.mul_monomial(&self.coefficients[0], &self.exponents);
}
if rhs.nterms() == 1 {
return self
.clone()
.mul_monomial(&rhs.coefficients[0], &rhs.exponents);
}
if let Some(r) = self.mul_dense(rhs) {
r
} else {
self.heap_mul(rhs)
}
}
}
impl<'a, F: Ring, E: Exponent> Mul<&'a MultivariatePolynomial<F, E, LexOrder>>
for MultivariatePolynomial<F, E, LexOrder>
{
type Output = MultivariatePolynomial<F, E, LexOrder>;
/// Multiply two polynomials, using either use dense multiplication or heap multiplication.
#[inline]
fn mul(self, rhs: &'a MultivariatePolynomial<F, E, LexOrder>) -> Self::Output {
(&self) * rhs
}
}
impl<'a, F: EuclideanDomain, E: PositiveExponent> Div<&'a MultivariatePolynomial<F, E, LexOrder>>
for &MultivariatePolynomial<F, E, LexOrder>
{
type Output = MultivariatePolynomial<F, E, LexOrder>;
fn div(self, other: &'a MultivariatePolynomial<F, E, LexOrder>) -> Self::Output {
self.try_div(other)
.unwrap_or_else(|| panic!("No clean division of {self} by {other}"))
}
}
impl<'a, F: EuclideanDomain, E: PositiveExponent> Div<&'a MultivariatePolynomial<F, E, LexOrder>>
for MultivariatePolynomial<F, E, LexOrder>
{
type Output = MultivariatePolynomial<F, E, LexOrder>;
fn div(
self: MultivariatePolynomial<F, E, LexOrder>,
other: &'a MultivariatePolynomial<F, E, LexOrder>,
) -> Self::Output {
(&self).div(other)
}
}
impl<F: Ring, E: Exponent, O: MonomialOrder> MultivariatePolynomial<F, E, O> {
/// Change the monomial order of the polynomial from `O` to `ON`.
pub fn reorder<ON: MonomialOrder>(&self) -> MultivariatePolynomial<F, E, ON> {
let mut sorted_index: Vec<_> = (0..self.nterms()).collect();
sorted_index.sort_by(|a, b| ON::cmp(self.exponents(*a), self.exponents(*b)));
let coefficients: Vec<_> = sorted_index
.iter()
.map(|i| self.coefficients[*i].clone())
.collect();
let exponents: Vec<_> = sorted_index
.iter()
.flat_map(|i| self.exponents(*i))
.cloned()
.collect();
MultivariatePolynomial {
coefficients,
exponents,
ring: self.ring.clone(),
variables: self.variables.clone(),
_phantom: PhantomData,
}
}
/// Multiply every coefficient with `other`.
pub fn mul_coeff(mut self, other: F::Element) -> Self {
if self.ring.is_one(&other) {
return self;
}
for c in &mut self.coefficients {
self.ring.mul_assign(c, &other);
}
for i in (0..self.nterms()).rev() {
if self.ring.is_zero(&self.coefficients[i]) {
self.coefficients.remove(i);
self.exponents
.drain(i * self.nvars()..(i + 1) * self.nvars());
}
}
self
}
/// Map a coefficient using the function `f`.
pub fn map_coeff<U: Ring, T: Fn(&F::Element) -> U::Element>(
&self,
f: T,
field: U,
) -> MultivariatePolynomial<U, E, O> {
let mut coefficients = Vec::with_capacity(self.coefficients.len());
let mut exponents = Vec::with_capacity(self.exponents.len());
for m in self.into_iter() {
let nc = f(m.coefficient);
if !field.is_zero(&nc) {
coefficients.push(nc);
exponents.extend(m.exponents);
}
}
MultivariatePolynomial {
coefficients,
exponents,
ring: field,
variables: self.variables.clone(),
_phantom: PhantomData,
}
}
/// Add `exponents` to every exponent.
pub fn mul_exp(mut self, exponents: &[E]) -> Self {
debug_assert_eq!(self.nvars(), exponents.len());
if self.nvars() == 0 {
return self;
}
for e in self.exponents_iter_mut() {
for (e1, e2) in e.iter_mut().zip(exponents) {
*e1 = e1.checked_add(e2).expect("overflow in adding exponents");
}
}
self
}
#[inline]
pub fn max_coeff(&self) -> &F::Element {
self.coefficients.last().unwrap()
}
#[inline]
pub fn max_exp(&self) -> &[E] {
if self.coefficients.is_empty() {
panic!("Cannot get max exponent of empty polynomial");
}
&self.exponents[(self.nterms() - 1) * self.nvars()..self.nterms() * self.nvars()]
}
/// Add a new monomial with coefficient `other` and exponent one.
pub fn add_constant(mut self, other: F::Element) -> Self {
let nvars = self.nvars();
self.append_monomial(other, &vec![E::zero(); nvars]);
self
}
#[inline]
pub fn mul_monomial(self, coefficient: &F::Element, exponents: &[E]) -> Self {
self.mul_coeff(coefficient.clone()).mul_exp(exponents)
}
/// Check if the polynomial contains a variable `x`.
pub fn contains(&self, x: usize) -> bool {
if self.nvars() == 0 {
return false;
}
for e in self.exponents.iter().skip(x).step_by(self.nvars()) {
if *e != E::zero() {
return true;
}
}
false
}
/// Get the degree of the variable `x`.
/// This operation is O(n).
pub fn degree(&self, x: usize) -> E {
if self.nvars() == 0 || self.is_zero() {
return E::zero();
}
let mut max = self.exponents[x];
for t in 1..self.nterms() {
let e = self.exponents[x + t * self.nvars()];
if max < e {
max = e;
}
}
max
}
/// Get the lowest and highest exponent of the variable `x`.
/// This operation is O(n).
pub fn degree_bounds(&self, x: usize) -> (E, E) {
if self.nvars() == 0 || self.is_zero() {
return (E::zero(), E::zero());
}
let mut min = None;
let mut max = None;
for e in self.exponents.iter().skip(x).step_by(self.nvars()) {
if max.map(|max| max < *e).unwrap_or(true) {
max = Some(*e);
}
if min.map(|min| min > *e).unwrap_or(true) {
min = Some(*e);
}
}
(min.unwrap_or(E::zero()), max.unwrap_or(E::zero()))
}
// Get the highest degree of a variable in the leading monomial.
pub fn ldegree(&self, v: usize) -> E {
if self.is_zero() {
return E::zero();
}
self.last_exponents()[v]
}
/// Get the highest degree of the leading monomial.
pub fn ldegree_max(&self) -> E {
if self.is_zero() {
return E::zero();
}
*self.last_exponents().iter().max().unwrap_or(&E::zero())
}
/// Get the leading coefficient.
pub fn lcoeff(&self) -> F::Element {
if self.is_zero() {
return self.ring.zero();
}
self.coefficients.last().unwrap().clone()
}
/// Perform self % var^pow.
pub fn mod_var(&self, var: usize, pow: E) -> Self {
let mut m = self.zero();
for t in self.into_iter() {
if t.exponents[var] < pow {
m.append_monomial(t.coefficient.clone(), t.exponents);
}
}
m
}
/// Take the derivative of the polynomial w.r.t the variable `var`.
pub fn derivative(&self, var: usize) -> Self {
debug_assert!(var < self.nvars());
let mut res = self.zero_with_capacity(self.nterms());
let mut exp = vec![E::zero(); self.nvars()];
for x in self {
if x.exponents[var] > E::zero() {
exp.copy_from_slice(x.exponents);
let pow = exp[var].to_i32() as u64;
exp[var] = exp[var] - E::one();
res.append_monomial(
self.ring.mul(x.coefficient, &self.ring.nth(pow.into())),
&exp,
);
}
}
res
}
/// Get the coefficient of the monomial with the given exponents if it is present.
pub fn coefficient(&self, exponents: &[E]) -> Option<F::Element> {
if self.is_zero() {
if exponents.iter().all(|e| *e == E::zero()) {
return Some(self.ring.zero());
}
return None;
}
let mut low = 0;
let mut high = self.coefficients.len();
while low < high {
let mid = low + (high - low) / 2;
match O::cmp(
&self.exponents[mid * self.nvars()..(mid + 1) * self.nvars()],
exponents,
) {
Ordering::Equal => return Some(self.coefficients[mid].clone()),
Ordering::Less => {
low = mid + 1;
}
Ordering::Greater => {
high = mid;
}
}
}
None
}
pub fn map_exp<E2: Exponent>(&self, f: impl Fn(&E) -> E2) -> MultivariatePolynomial<F, E2, O> {
MultivariatePolynomial {
coefficients: self.coefficients.clone(),
exponents: self.exponents.iter().map(f).collect::<Vec<_>>(),
ring: self.ring.clone(),
variables: self.variables.clone(),
_phantom: PhantomData,
}
}
/// Kronecker map all variables starting from `start_index` using the powers given in `powers`:
/// `x_i -> x_{start_index}^powers[i - start_index]` for `i > start_index`.
pub fn kronecker_map(&self, powers: &[E], start_index: usize) -> Self {
let mut res = self.zero_with_capacity(self.nterms());
let mut new_exponents = vec![E::zero(); self.nvars()];
for a in self {
new_exponents[..(start_index + 1)].copy_from_slice(&a.exponents[..(start_index + 1)]);
for (i, e) in a.exponents.iter().skip(start_index + 1).enumerate() {
new_exponents[start_index] += powers[i] * *e;
}
res.append_monomial(a.coefficient.clone(), &new_exponents);
}
res
}
/// Invert a Kronecker map.
pub fn kronecker_inv_map(&self, powers: &[E], start_index: usize) -> Self {
let mut res = self.zero_with_capacity(self.nterms());
let mut new_exponents = vec![E::zero(); self.nvars()];
for a in self {
new_exponents[..start_index].copy_from_slice(&a.exponents[..start_index]);
let mut total = a.exponents[start_index];
new_exponents[start_index] = total % powers[0];
for i in (start_index + 1)..self.nvars() {
let previous_power = if i == start_index + 1 {
E::one()
} else {
powers[i - start_index - 2]
};
total = total - new_exponents[i - 1] * previous_power;
new_exponents[i] = total % powers[i - start_index];
new_exponents[i] = new_exponents[i] / powers[i - start_index - 1];
}
res.append_monomial(a.coefficient.clone(), &new_exponents);
}
res
}
/// Create a polynomial from an unordered list of coefficients and flattened exponents.
pub fn from_coefficient_list(
mut coefficients: Vec<F::Element>,
exponents: Vec<E>,
vars: Arc<Vec<PolyVariable>>,
ring: &F,
) -> Self {
let nterms = coefficients.len();
let nvars = exponents.len() / coefficients.len();
let mut indices = (0..nterms).collect::<Vec<_>>();
indices.sort_unstable_by(|&i, &j| {
O::cmp(
&exponents[i * nvars..(i + 1) * nvars],
&exponents[j * nvars..(j + 1) * nvars],
)
});
let mut poly = MultivariatePolynomial::new(ring, Some(nterms), vars);
for i in indices {
poly.append_monomial_back(
std::mem::replace(&mut coefficients[i], ring.zero()),
&exponents[i * nvars..(i + 1) * nvars],
);
}
poly
}
}
impl<F: Ring, E: PositiveExponent> MultivariatePolynomial<F, E, LexOrder> {
/// Remove all non-occurring variables from the polynomial.
pub fn condense(&mut self) {
if self.nvars() == 0 {
return;
}
let degrees: Vec<_> = (0..self.nvars())
.filter(|i| self.degree(*i) > E::zero())
.collect();
let mut new_exponents = vec![E::zero(); self.nterms() * degrees.len()];
if degrees.is_empty() {
self.exponents = new_exponents;
self.variables = Arc::new(vec![]);
return;
}
for (d, e) in new_exponents
.chunks_mut(degrees.len())
.zip(self.exponents_iter())
{
for (dr, s) in d.iter_mut().zip(°rees) {
*dr = e[*s];
}
}
self.exponents = new_exponents;
self.variables = Arc::new(
degrees
.into_iter()
.map(|x| self.variables[x].clone())
.collect(),
);
}
/// Replace a variable `n` in the polynomial by an element from
/// the ring `v`.
pub fn replace(&self, n: usize, v: &F::Element) -> MultivariatePolynomial<F, E, LexOrder> {
if (n + 1..self.nvars()).all(|i| self.degree(i) == E::zero()) {
return self.replace_last(n, v);
}
let mut coefficients = self.coefficients.clone();
let mut exponents = self.exponents.clone();
// TODO: cache power taking?
for (coefficient, exponent) in coefficients
.iter_mut()
.zip(exponents.chunks_mut(self.nvars()))
{
if exponent[n] == E::zero() {
continue;
}
self.ring
.mul_assign(coefficient, &self.ring.pow(v, exponent[n].to_i32() as u64));
exponent[n] = E::zero();
}
Self::from_coefficient_list(coefficients, exponents, self.variables.clone(), &self.ring)
}
/// Replace the last variable `n` in the polynomial by an element from
/// the ring `v`.
pub fn replace_last(&self, n: usize, v: &F::Element) -> MultivariatePolynomial<F, E, LexOrder> {
const MAX_EXP_BUF: usize = 100000;
let mut res = self.zero_with_capacity(self.nterms());
let mut e = vec![E::zero(); self.nvars()];
let mut res_cache =
vec![self.ring.zero(); (self.degree(n).to_i32() as usize + 1).min(MAX_EXP_BUF)];
let nvars = self.nvars();
for t in self {
if t.exponents[n] == E::zero() {
res.append_monomial(t.coefficient.clone(), t.exponents);
continue;
}
let ee = t.exponents[n].to_i32() as usize;
let c = if ee < MAX_EXP_BUF {
if self.ring.is_zero(&res_cache[ee]) {
res_cache[ee] = self.ring.pow(v, ee as u64);
}
self.ring.mul(t.coefficient, &res_cache[ee])
} else {
self.ring.mul(t.coefficient, &self.ring.pow(v, ee as u64))
};
if self.ring.is_zero(&c) {
continue;
}
e.copy_from_slice(t.exponents);
e[n] = E::zero();
let mut different = false;
if !res.is_zero() {
let start = (res.nterms() - 1) * nvars;
for i in (0..n).rev() {
if unsafe { e.get_unchecked(i) != res.exponents.get_unchecked(start + i) } {
different = true;
break;
}
}
} else {
different = true;
}
if different {
res.coefficients.push(c);
res.exponents.extend_from_slice(&e);
} else {
let l = res.coefficients.last_mut().unwrap();
self.ring.add_assign(l, &c);
if self.ring.is_zero(l) {
res.coefficients.pop();
res.exponents.truncate(res.exponents.len() - self.nvars());
}
}
}
res
}
pub fn evaluate<T: FloatLike, M: Fn(&F::Element) -> T>(&self, map_coeff: M, point: &[T]) -> T {
let mut res = map_coeff(&self.ring.zero());
for t in self {
let mut c = map_coeff(t.coefficient);
for (i, v) in point.iter().zip(t.exponents) {
if v != &E::zero() {
c *= i.pow(v.to_u32() as u64);
}
}
res += c;
}
res
}
/// Evaluate the polynomial at the given point, mapping coefficients to the ring `U`.
pub fn evaluate_with_coeff_map<U: Ring, T: Fn(&F::Element) -> U::Element>(
&self,
map_coeff: T,
point: &[U::Element],
ring: &U,
) -> U::Element {
let mut res = map_coeff(&self.ring.zero());
assert_eq!(point.len(), self.nvars());
for t in self {
let mut c = map_coeff(t.coefficient);
for (i, v) in point.iter().zip(t.exponents) {
if v != &E::zero() {
ring.mul_assign(&mut c, &ring.pow(i, v.to_u32() as u64));
}
}
ring.add_assign(&mut res, &c);
}
res
}
/// Replace all variables in the polynomial by an element from
/// the ring `v`.
pub fn replace_all(&self, r: &[F::Element]) -> F::Element {
let mut res = self.ring.zero();
// TODO: cache power taking?
for t in self {
let mut c = t.coefficient.clone();
for (i, v) in r.iter().zip(t.exponents) {
if v != &E::zero() {
self.ring
.mul_assign(&mut c, &self.ring.pow(i, v.to_u32() as u64));
}
}
self.ring.add_assign(&mut res, &c);
}
res
}
/// Replace a variable `n` in the polynomial by a polynomial `v`.
pub fn replace_with_poly(&self, n: usize, v: &Self) -> Self {
assert_eq!(self.variables, v.variables);
if v.is_constant() {
return self.replace(n, &v.lcoeff());
}
let mut res = self.zero_with_capacity(self.nterms());
let mut exp = vec![E::zero(); self.nvars()];
for t in self {
if t.exponents[n] == E::zero() {
res.append_monomial(t.coefficient.clone(), &t.exponents[..self.nvars()]);
continue;
}
exp.copy_from_slice(t.exponents);
exp[n] = E::zero();
// TODO: cache v^e
res = res
+ (&v.pow(t.exponents[n].to_i32() as usize)
* &self.monomial(t.coefficient.clone(), exp.clone()));
}
res
}
/// Replace all variables except `v` in the polynomial by elements from
/// the ring.
pub fn replace_except(
&self,
v: usize,
r: &[(usize, F::Element)],
cache: &mut [Vec<F::Element>],
) -> MultivariatePolynomial<F, E, LexOrder> {
let mut tm: HashMap<E, F::Element> = HashMap::new();
for t in self {
let mut c = t.coefficient.clone();
for (n, vv) in r {
let p = t.exponents[*n].to_i32() as usize;
if p > 0 {
if p < cache[*n].len() {
if self.ring.is_zero(&cache[*n][p]) {
cache[*n][p] = self.ring.pow(vv, p as u64);
}
self.ring.mul_assign(&mut c, &cache[*n][p]);
} else {
self.ring.mul_assign(&mut c, &self.ring.pow(vv, p as u64));
}
}
}
tm.entry(t.exponents[v])
.and_modify(|e| self.ring.add_assign(e, &c))
.or_insert(c);
}
let mut res = self.zero();
let mut e = vec![E::zero(); self.nvars()];
for (k, c) in tm {
e[v] = k;
res.append_monomial(c, &e);
e[v] = E::zero();
}
res
}
/// Shift a variable `var` to `var+shift`.
pub fn shift_var(&self, var: usize, shift: &F::Element) -> Self {
let d = self.degree(var).to_i32() as usize;
let y_poly = self.to_univariate_polynomial_list(var);
let mut v = vec![self.zero(); d + 1];
for (x_poly, p) in y_poly {
v[p.to_i32() as usize] = x_poly;
}
for k in 0..d {
for j in (k..d).rev() {
v[j] = &v[j] + &v[j + 1].clone().mul_coeff(shift.clone());
}
}
let mut shifted_coefficients = Vec::with_capacity(v.len());
for (i, mut v) in v.into_iter().enumerate() {
for x in v.exponents.chunks_mut(self.nvars()) {
x[var] = E::from_i32(i as i32);
}
shifted_coefficients.push(v);
}
self.merge_shifted_univariate_coefficients(shifted_coefficients)
}
fn merge_shifted_univariate_coefficients(&self, coefficients: Vec<Self>) -> Self {
let capacity = coefficients.iter().map(|p| p.nterms()).sum();
let mut poly = self.zero_with_capacity(capacity);
let mut indices = vec![0usize; coefficients.len()];
let mut exponent: SmallVec<[E; INLINED_EXPONENTS]> = smallvec![E::zero(); self.nvars()];
loop {
let mut min_poly: Option<usize> = None;
for (i, p) in coefficients.iter().enumerate() {
if indices[i] == p.nterms() {
continue;
}
if min_poly
.map(|j| {
LexOrder::cmp(
p.exponents(indices[i]),
coefficients[j].exponents(indices[j]),
)
.is_lt()
})
.unwrap_or(true)
{
min_poly = Some(i);
}
}
let Some(min_poly) = min_poly else {
return poly;
};
exponent.clear();
exponent.extend_from_slice(coefficients[min_poly].exponents(indices[min_poly]));
let mut coefficient = self.ring.zero();
for (i, p) in coefficients.iter().enumerate() {
if indices[i] < p.nterms() && p.exponents(indices[i]) == &exponent[..] {
self.ring
.add_assign(&mut coefficient, p.coefficients[indices[i]].clone());
indices[i] += 1;
}
}
poly.append_monomial_back(coefficient, &exponent);
}
}
/// Compute the inverse of the univariate polynomial in `var` up until `pow` using Newton's method.
pub fn inverse_univariate(&self, var: usize, pow: E) -> Self {
let mut g = self.constant(
self.ring
.try_div(&self.ring.one(), &self.get_constant())
.unwrap(),
);
let mut exp = E::one();
while exp < pow {
exp = exp * E::from_u32(2);
let h =
g.clone().mul_coeff(self.ring.nth(2.into())) - self * &(&g * &g).mod_var(var, exp);
g = h.mod_var(var, exp);
}
g
}
}
impl<F: Ring, E: Exponent> MultivariatePolynomial<F, E, LexOrder> {
/// Check if all exponents are positive.
pub fn is_polynomial(&self) -> bool {
self.is_zero() || self.exponents.iter().all(|e| *e >= E::zero())
}
/// Get the leading coefficient under a given variable ordering.
/// This operation is O(n) if the variables are out of order.
pub fn lcoeff_varorder(&self, vars: &[usize]) -> F::Element {
if vars.windows(2).all(|s| s[0] < s[1]) {
return self.lcoeff();
}
let mut highest = vec![E::zero(); self.nvars()];
let mut highestc = &self.ring.zero();
'nextmon: for m in self.into_iter() {
let mut more = false;
for &v in vars {
if more {
highest[v] = m.exponents[v];
} else {
match m.exponents[v].cmp(&highest[v]) {
Ordering::Less => {
continue 'nextmon;
}
Ordering::Greater => {
highest[v] = m.exponents[v];
more = true;
}
Ordering::Equal => {}
}
}
}
highestc = m.coefficient;
}
debug_assert!(!self.ring.is_zero(highestc));
highestc.clone()
}
/// Get the leading coefficient of a multivariate polynomial viewed as a
/// univariate polynomial in `x`.
pub fn univariate_lcoeff(&self, x: usize) -> MultivariatePolynomial<F, E, LexOrder> {
let d = self.degree(x);
let mut lcoeff = self.zero();
if self.coefficients.is_empty() {
return lcoeff;
}
if d == E::zero() {
return self.clone();
}
let mut e = vec![E::zero(); self.nvars()];
for t in self {
if t.exponents[x] == d {
e.copy_from_slice(t.exponents);
e[x] = E::zero();
lcoeff.append_monomial(t.coefficient.clone(), &e);
}
}
lcoeff
}
/// Count the number of terms that have the maximum degree in the given variable.
pub fn terms_with_max_degree(&self, x: usize) -> usize {
let mut max_degree = self.exponents[x];
let mut count = 0;
for t in 1..self.nterms() {
let e = self.exponents[t * self.nvars() + x];
if e == max_degree {
count += 1;
} else if e > max_degree {
max_degree = e;
count = 1;
}
}
count
}
/// Get the bivariate degree of a multivariate polynomial viewed as a bivariate polynomial in the first two variables.
pub fn bivariate_deg(&self) -> (E, E) {
if self.is_zero() {
return (E::zero(), E::zero());
}
let mut d = self.exponents[0];
let mut d1 = self.exponents[1];
let nvars = self.nvars();
assert!(self.exponents.len() == self.nterms() * nvars);
for e in 0..self.nterms() {
if self.exponents[e * nvars] > d {
d = self.exponents[e * nvars];
d1 = self.exponents[e * nvars + 1];
} else if self.exponents[e * nvars] == d && self.exponents[e * nvars + 1] > d1 {
d1 = self.exponents[e * nvars + 1];
}
}
(d, d1)
}
/// Get the leading coefficient of a multivariate polynomial viewed as a bivariate polynomial in the first two variables.
pub fn bivariate_lcoeff(&self) -> MultivariatePolynomial<F, E, LexOrder> {
let mut lcoeff = self.zero();
if self.coefficients.is_empty() {
return lcoeff;
}
let (d, d1) = self.bivariate_deg();
let mut e = vec![E::zero(); self.nvars()];
for t in self {
if t.exponents[0] == d && t.exponents[1] == d1 {
e.copy_from_slice(t.exponents);
e[0] = E::zero();
e[1] = E::zero();
lcoeff.append_monomial(t.coefficient.clone(), &e);
}
}
lcoeff
}
/// Get the leading coefficient viewed as a polynomial
/// in all variables except the last variable `n`.
pub fn lcoeff_last(&self, n: usize) -> MultivariatePolynomial<F, E, LexOrder> {
if self.is_zero() {
return self.clone();
}
// the last variable should have the least sorting priority,
// so the last term should still be the lcoeff
let last = self.last_exponents();
let mut res = self.zero();
let mut e: SmallVec<[E; INLINED_EXPONENTS]> = smallvec![E::zero(); self.nvars()];
for t in (0..self.nterms()).rev() {
if (0..self.nvars() - 1).all(|i| self.exponents(t)[i] == last[i] || i == n) {
e[n] = self.exponents(t)[n];
res.append_monomial(self.coefficients[t].clone(), &e);
e[n] = E::zero();
} else {
break;
}
}
res
}
/// Get the leading coefficient viewed as a polynomial
/// in all variables with order as described in `vars` except the last variable in `vars`.
/// This operation is O(n) if the variables are out of order.
pub fn lcoeff_last_varorder(&self, vars: &[usize]) -> MultivariatePolynomial<F, E, LexOrder> {
if self.is_zero() {
return self.clone();
}
if vars.windows(2).all(|s| s[0] < s[1]) {
return self.lcoeff_last(*vars.last().unwrap());
}
let (vars, lastvar) = vars.split_at(vars.len() - 1);
let mut highest = vec![E::zero(); self.nvars()];
let mut indices = Vec::with_capacity(10);
'nextmon: for (i, m) in self.into_iter().enumerate() {
let mut more = false;
for &v in vars {
if more {
highest[v] = m.exponents[v];
} else {
match m.exponents[v].cmp(&highest[v]) {
Ordering::Less => {
continue 'nextmon;
}
Ordering::Greater => {
highest[v] = m.exponents[v];
indices.clear();
more = true;
}
Ordering::Equal => {}
}
}
}
indices.push(i);
}
let mut res = self.zero();
let mut e = vec![E::zero(); self.nvars()];
for i in indices {
e[lastvar[0]] = self.exponents(i)[lastvar[0]];
res.append_monomial(self.coefficients[i].clone(), &e);
e[lastvar[0]] = E::zero();
}
res
}
/// Change the order of the variables in the polynomial, using `order`.
/// The map can also be reversed, by setting `inverse` to `true`.
pub(crate) fn rearrange_impl(
&self,
order: &[usize],
inverse: bool,
update_variables: bool,
) -> MultivariatePolynomial<F, E, LexOrder> {
let mut new_exp = vec![E::zero(); self.nterms() * self.nvars()];
for (e, er) in new_exp.chunks_mut(self.nvars()).zip(self.exponents_iter()) {
for x in 0..order.len() {
if !inverse {
e[x] = er[order[x]];
} else {
e[order[x]] = er[x];
}
}
}
let mut indices: Vec<usize> = (0..self.nterms()).collect();
indices.sort_unstable_by_key(|&i| &new_exp[i * self.nvars()..(i + 1) * self.nvars()]);
let mut res = self.zero_with_capacity(self.nterms());
for i in indices {
res.append_monomial(
self.coefficients[i].clone(),
&new_exp[i * self.nvars()..(i + 1) * self.nvars()],
);
}
if update_variables {
let mut vm = self.variables.as_ref().clone();
for x in 0..order.len() {
if !inverse {
vm[x] = self.variables[order[x]].clone();
} else {
vm[order[x]] = self.variables[x].clone();
}
}
res.variables = Arc::new(vm);
}
res
}
/// Change the order of the variables in the polynomial, using `order`.
/// The map can also be reversed, by setting `inverse` to `true`.
pub fn rearrange(
&self,
order: &[usize],
inverse: bool,
) -> MultivariatePolynomial<F, E, LexOrder> {
self.rearrange_impl(order, inverse, true)
}
/// Change the order of the variables in the polynomial, using `order`.
/// The order may contain `None`, to signal unmapped indices. This operation
/// allows the polynomial to grow in size.
pub fn rearrange_with_growth(
&self,
order: &[PolyVariable],
) -> Result<MultivariatePolynomial<F, E, LexOrder>, String> {
let new_order: Vec<_> = order
.iter()
.map(|x| self.variables.iter().position(|e| e == x))
.collect();
for (i, v) in self.variables.iter().enumerate() {
if !new_order.contains(&Some(i)) && self.contains(i) {
return Err(format!("Variable {v} is not in the new order"));
}
}
let mut new_exp = vec![E::zero(); self.nterms() * order.len()];
for (e, er) in new_exp.chunks_mut(order.len()).zip(self.exponents_iter()) {
for x in 0..new_order.len() {
if let Some(v) = new_order[x] {
e[x] = er[v];
}
}
}
let mut indices: Vec<usize> = (0..self.nterms()).collect();
indices.sort_unstable_by_key(|&i| &new_exp[i * order.len()..(i + 1) * order.len()]);
let mut res =
MultivariatePolynomial::new(&self.ring, self.nterms().into(), Arc::new(order.to_vec()));
for i in indices {
res.append_monomial(
self.coefficients[i].clone(),
&new_exp[i * order.len()..(i + 1) * order.len()],
);
}
Ok(res)
}
/// Compute `self^pow`.
pub fn pow(&self, pow: usize) -> Self {
if pow == 0 {
return self.one();
}
if pow == 1 {
return self.clone();
}
if self.is_constant() {
return self.constant(self.ring.pow(&self.lcoeff(), pow as u64));
}
if self.coefficients.len() == 1 {
return self.monomial(
self.ring.pow(&self.coefficients[0], pow as u64),
self.exponents
.iter()
.map(|x| *x * E::from_i32(pow as i32))
.collect(),
);
}
// heuristic for when to use heap_pow
if pow > 10 || (0..self.nvars()).all(|x| self.degree(x) <= E::from_i32(2)) {
// if the characteristic is non-zero, a division by the exponent in the heap_pow algorithm
// may cause a division by 0
if self.ring.characteristic() == 0
|| self.nvars() == 1
&& self.degree(0).to_i32() as usize + 1 < self.ring.characteristic()
{
return self.heap_pow(pow);
}
}
// perform repeated multiplication instead of binary exponentiation, as
// the latter is often much slower for sparse polynomials
let mut res = self * self;
for _ in 2..pow {
res = &res * self;
}
res
}
pub fn to_univariate(&self, var: usize) -> UnivariatePolynomial<PolynomialRing<F, E>> {
let c = self.to_univariate_polynomial_list(var);
let mut p = UnivariatePolynomial::new(
&PolynomialRing::from_poly(self),
None,
Arc::new(self.variables[var].clone()),
);
if c.is_empty() {
return p;
}
p.coefficients = vec![self.zero(); c.last().unwrap().1.to_i32() as usize + 1];
for (q, e) in c {
if e < E::zero() {
panic!("Negative exponent in univariate conversion");
}
p.coefficients[e.to_i32() as usize] = q;
}
p
}
pub fn to_univariate_from_univariate(&self, var: usize) -> UnivariatePolynomial<F> {
let mut p =
UnivariatePolynomial::new(&self.ring, None, Arc::new(self.variables[var].clone()));
if self.is_zero() {
return p;
}
p.coefficients = vec![p.ring.zero(); self.degree(var).to_i32() as usize + 1];
for (q, e) in self.coefficients.iter().zip(self.exponents_iter()) {
if e[var] < E::zero() {
panic!("Negative exponent in univariate conversion");
}
p.coefficients[e[var].to_i32() as usize] = q.clone();
}
p
}
/// Create a univariate polynomial coefficient list out of a multivariate polynomial.
/// The output is sorted in the degree.
pub fn to_univariate_polynomial_list(
&self,
x: usize,
) -> Vec<(MultivariatePolynomial<F, E, LexOrder>, E)> {
if self.coefficients.is_empty() {
return vec![];
}
// get maximum degree for variable x
let mut mindeg = E::zero();
let mut maxdeg = E::zero();
for t in 0..self.nterms() {
let d = self.exponents(t)[x];
if d > maxdeg {
maxdeg = d;
}
if d < mindeg {
mindeg = d;
}
}
// construct the coefficient per power of x
let mut result = vec![];
let mut e: SmallVec<[E; INLINED_EXPONENTS]> = smallvec![E::zero(); self.nvars()];
for d in mindeg.to_i32()..maxdeg.to_i32() + 1 {
// TODO: add bounds estimate
let mut a = self.zero();
for t in 0..self.nterms() {
if self.exponents(t)[x].to_i32() == d {
for (i, ee) in self.exponents(t).iter().enumerate() {
e[i] = *ee;
}
e[x] = E::zero();
a.append_monomial(self.coefficients[t].clone(), &e);
}
}
if !a.is_zero() {
result.push((a, E::from_i32(d)));
}
}
result
}
/// Split the polynomial as a polynomial in `xs` if include is true,
/// else excluding `xs`.
pub fn to_multivariate_polynomial_list(
&self,
xs: &[usize],
include: bool,
) -> HashMap<SmallVec<[E; INLINED_EXPONENTS]>, MultivariatePolynomial<F, E, LexOrder>> {
if self.coefficients.is_empty() {
return HashMap::new();
}
let mut tm: HashMap<
SmallVec<[E; INLINED_EXPONENTS]>,
MultivariatePolynomial<F, E, LexOrder>,
> = HashMap::new();
let mut e_not_in_xs = smallvec![E::zero(); self.nvars()];
let mut e_in_xs = smallvec![E::zero(); self.nvars()];
for t in self {
for (i, ee) in t.exponents.iter().enumerate() {
e_not_in_xs[i] = *ee;
e_in_xs[i] = E::zero();
}
for x in xs {
e_in_xs[*x] = e_not_in_xs[*x];
e_not_in_xs[*x] = E::zero();
}
if include {
tm.entry(e_in_xs.clone())
.and_modify(|x| x.append_monomial(t.coefficient.clone(), &e_not_in_xs))
.or_insert_with(|| {
MultivariatePolynomial::monomial(
self,
t.coefficient.clone(),
e_not_in_xs.to_vec(),
)
});
} else {
tm.entry(e_not_in_xs.clone())
.and_modify(|x| x.append_monomial(t.coefficient.clone(), &e_in_xs))
.or_insert_with(|| {
MultivariatePolynomial::monomial(
self,
t.coefficient.clone(),
e_in_xs.to_vec(),
)
});
}
}
tm
}
pub(crate) fn mul_univariate_dense(&self, rhs: &Self, max_pow: Option<usize>) -> Self {
if self.is_constant() {
if let Some(m) = max_pow
&& let Some(var) = rhs.last_exponents().iter().position(|e| *e != E::zero())
&& rhs.degree(var).to_i32() > m as i32
{
return rhs
.mod_var(var, E::from_i32(m as i32 + 1))
.mul_coeff(self.lcoeff());
}
return rhs.clone().mul_coeff(self.lcoeff());
}
if rhs.is_constant() {
if let Some(m) = max_pow
&& let Some(var) = self.last_exponents().iter().position(|e| *e != E::zero())
&& self.degree(var).to_i32() > m as i32
{
return self
.mod_var(var, E::from_i32(m as i32 + 1))
.mul_coeff(rhs.lcoeff());
}
return self.clone().mul_coeff(rhs.lcoeff());
}
let var = self
.last_exponents()
.iter()
.position(|e| *e != E::zero())
.unwrap();
let d1 = self.degree(var);
let d2 = rhs.degree(var);
let mut max = (d1.to_i32() + d2.to_i32()) as usize;
if let Some(m) = max_pow {
max = max.min(m);
}
let mut coeffs = vec![self.ring.zero(); max + 1];
for x in self {
for y in rhs {
let pos = x.exponents[var].to_i32() + y.exponents[var].to_i32();
if pos as usize > max {
continue;
}
self.ring
.add_mul_assign(&mut coeffs[pos as usize], x.coefficient, y.coefficient);
}
}
let mut exp = vec![E::zero(); self.nvars()];
let mut res = self.zero_with_capacity(coeffs.len());
for (p, c) in coeffs.into_iter().enumerate() {
if !self.ring.is_zero(&c) {
exp[var] = E::from_i32(p as i32);
res.append_monomial(c, &exp);
}
}
res
}
/// Synthetic division for univariate polynomials, where `div` is monic.
pub(crate) fn quot_rem_univariate_monic(
&self,
div: &MultivariatePolynomial<F, E, LexOrder>,
) -> (
MultivariatePolynomial<F, E, LexOrder>,
MultivariatePolynomial<F, E, LexOrder>,
) {
debug_assert_eq!(div.lcoeff(), self.ring.one());
if self.is_zero() {
return (self.clone(), self.clone());
}
let mut dividendpos = self.nterms() - 1; // work from the back
let mut q = self.zero_with_capacity(self.nterms());
let mut r = self.zero();
// determine the variable
let mut var = 0;
for (i, x) in self.last_exponents().iter().enumerate() {
if !x.is_zero() {
var = i;
break;
}
}
let m = div.ldegree_max();
let mut pow = self.ldegree_max();
loop {
// find the power in the dividend if it exists
let mut coeff = loop {
if self.exponents(dividendpos)[var] == pow {
break self.coefficients[dividendpos].clone();
}
if dividendpos == 0 || self.exponents(dividendpos)[var] < pow {
break self.ring.zero();
}
dividendpos -= 1;
};
let mut qindex = 0; // starting from highest
let mut bindex = 0; // starting from lowest
while bindex < div.nterms() && qindex < q.nterms() {
while bindex + 1 < div.nterms()
&& div.exponents(bindex)[var] + q.exponents(qindex)[var] < pow
{
bindex += 1;
}
if div.exponents(bindex)[var] + q.exponents(qindex)[var] == pow {
self.ring.sub_mul_assign(
&mut coeff,
&div.coefficients[bindex],
&q.coefficients[qindex],
);
}
qindex += 1;
}
if !self.ring.is_zero(&coeff) {
// can the division be performed? if not, add to rest
// TODO: refactor
let (quot, div) = if pow >= m {
(coeff, true)
} else {
(coeff, false)
};
if div {
let nterms = q.nterms();
let nvars = q.nvars();
q.coefficients.push(quot);
q.exponents.resize((nterms + 1) * nvars, E::zero());
q.exponents[nterms * nvars + var] = pow - m;
} else {
let nterms = r.nterms();
let nvars = r.nvars();
r.coefficients.push(quot);
r.exponents.resize((nterms + 1) * nvars, E::zero());
r.exponents[nterms * nvars + var] = pow;
}
}
if pow.is_zero() {
break;
}
pow = pow - E::one();
}
q.reverse_monomials();
r.reverse_monomials();
#[cfg(debug_assertions)]
{
if !(&q * div + r.clone() - self.clone()).is_zero() {
panic!("Division failed: ({self})/({div}): q={q}, r={r}");
}
}
(q, r)
}
fn mul_dense(
&self,
rhs: &MultivariatePolynomial<F, E, LexOrder>,
) -> Option<MultivariatePolynomial<F, E, LexOrder>> {
if !self.is_polynomial() || !rhs.is_polynomial() {
return None;
}
let max_degs_rev = (0..self.nvars())
.rev()
.map(|i| 1 + self.degree(i).to_i32() as usize + rhs.degree(i).to_i32() as usize)
.collect::<Vec<_>>();
if max_degs_rev.iter().filter(|x| **x > 1).count() == 1 {
if max_degs_rev.iter().sum::<usize>() < 10000 {
return Some(self.mul_univariate_dense(rhs, None));
}
return None;
}
let mut total: usize = 1;
for x in &max_degs_rev {
if *x > MAX_DENSE_MUL_BUFFER_SIZE {
return None;
}
if let Some(r) = total.checked_mul(*x) {
total = r;
} else {
return None;
}
}
if total > MAX_DENSE_MUL_BUFFER_SIZE {
return None;
}
#[inline(always)]
fn to_uni_var<E: Exponent>(s: &[E], max_degs_rev: &[usize]) -> u32 {
let mut shift = 1;
let mut res = s.last().unwrap().to_i32() as u32;
for (ee, &x) in s.iter().rev().skip(1).zip(max_degs_rev) {
shift *= x as u32;
res += ee.to_i32() as u32 * shift;
}
res
}
#[inline(always)]
fn from_uni_var<E: Exponent>(mut p: u32, max_degs_rev: &[usize], exp: &mut [E]) {
for (ee, &x) in exp.iter_mut().rev().zip(max_degs_rev) {
*ee = E::from_i32((p % x as u32) as i32);
p /= x as u32;
}
}
let mut uni_exp_self = vec![0; self.coefficients.len()];
for (es, s) in &mut uni_exp_self.iter_mut().zip(self.exponents_iter()) {
*es = to_uni_var(s, &max_degs_rev);
}
let mut uni_exp_rhs = vec![0; rhs.coefficients.len()];
for (es, s) in &mut uni_exp_rhs.iter_mut().zip(rhs.exponents_iter()) {
*es = to_uni_var(s, &max_degs_rev);
}
let mut exp = vec![E::zero(); self.nvars()];
let mut r = self.zero_with_capacity(self.nterms().max(rhs.nterms()));
// check if we need to use a dense indexing array to save memory
if total < 1000 {
let mut coeffs = vec![self.ring.zero(); total];
for (c1, e1) in self.coefficients.iter().zip(&uni_exp_self) {
for (c2, e2) in rhs.coefficients.iter().zip(&uni_exp_rhs) {
let pos = *e1 as usize + *e2 as usize;
self.ring.add_mul_assign(&mut coeffs[pos], c1, c2);
}
}
for (p, c) in coeffs.into_iter().enumerate() {
if !self.ring.is_zero(&c) {
from_uni_var(p as u32, &max_degs_rev, &mut exp);
r.append_monomial(c, &exp);
}
}
Some(r)
} else {
let mut coeffs = Vec::with_capacity(self.nterms().max(rhs.nterms()));
let mut coeff_index = DENSE_MUL_BUFFER.take();
if coeff_index.len() < total {
coeff_index.resize(total, 0u32);
}
for (c1, e1) in self.coefficients.iter().zip(&uni_exp_self) {
for (c2, e2) in rhs.coefficients.iter().zip(&uni_exp_rhs) {
let pos = *e1 as usize + *e2 as usize;
if coeff_index[pos] == 0 {
coeffs.push(self.ring.mul(c1, c2));
coeff_index[pos] = coeffs.len() as u32;
} else {
self.ring.add_mul_assign(
&mut coeffs[coeff_index[pos] as usize - 1],
c1,
c2,
);
}
}
}
for (p, c) in coeff_index[..total].iter_mut().enumerate() {
if *c != 0 {
from_uni_var(p as u32, &max_degs_rev, &mut exp);
r.append_monomial(
std::mem::replace(&mut coeffs[*c as usize - 1], self.ring.zero()),
&exp,
);
*c = 0;
}
}
DENSE_MUL_BUFFER.set(coeff_index);
Some(r)
}
}
fn heap_mul(
&self,
rhs: &MultivariatePolynomial<F, E, LexOrder>,
) -> MultivariatePolynomial<F, E, LexOrder> {
// place the smallest polynomial first, as this is faster
// in the heap algorithm
if self.nterms() > rhs.nterms() {
return rhs.heap_mul(self);
}
let degree_sum: Vec<_> = (0..self.nvars())
.map(|i| self.degree(i).to_i32() as i64 + rhs.degree(i).to_i32() as i64)
.collect();
// use a special routine if the exponents can be packed into a u64
let mut pack_u8 = true;
if self.nvars() <= 8
&& self.is_polynomial()
&& rhs.is_polynomial()
&& degree_sum.iter().all(|deg| {
if *deg > 255 {
pack_u8 = false;
}
*deg <= 255 || self.nvars() <= 4 && *deg <= 65535
})
{
return self.heap_mul_packed_exp(rhs, pack_u8);
}
let mut monomials = Vec::with_capacity(self.nterms() * self.nvars());
monomials.extend(
self.exponents(0)
.iter()
.zip(rhs.exponents(0))
.map(|(e1, e2)| *e1 + *e2),
);
let monomials = UnsafeCell::new((self.nvars(), monomials));
/// In order to prevent allocations of the exponents, store them in a single
/// append-only vector and use a key to index into it. For performance,
/// we use an unsafe cell.
#[derive(Clone, Copy)]
struct Key<'a, E: Exponent> {
index: usize,
monomials: &'a UnsafeCell<(usize, Vec<E>)>,
}
impl<E: Exponent> PartialEq for Key<'_, E> {
#[inline(always)]
fn eq(&self, other: &Self) -> bool {
unsafe {
let b1 = &*self.monomials.get();
b1.1.get_unchecked(self.index..self.index + b1.0)
== b1.1.get_unchecked(other.index..other.index + b1.0)
}
}
}
impl<E: Exponent> Eq for Key<'_, E> {}
impl<E: Exponent> PartialOrd for Key<'_, E> {
#[inline(always)]
fn partial_cmp(&self, other: &Self) -> Option<Ordering> {
Some(self.cmp(other))
}
}
impl<E: Exponent> Ord for Key<'_, E> {
#[inline(always)]
fn cmp(&self, other: &Self) -> Ordering {
unsafe {
let b1 = &*self.monomials.get();
b1.1.get_unchecked(self.index..self.index + b1.0)
.cmp(b1.1.get_unchecked(other.index..other.index + b1.0))
}
}
}
impl<E: Exponent> std::hash::Hash for Key<'_, E> {
#[inline(always)]
fn hash<H: std::hash::Hasher>(&self, state: &mut H) {
unsafe {
let b = &*self.monomials.get();
b.1.get_unchecked(self.index..self.index + b.0).hash(state);
}
}
}
let mut res = self.zero_with_capacity(self.nterms().max(rhs.nterms()));
let mut cache: HashMap<_, Vec<(usize, usize)>> = HashMap::new();
let mut q_cache: Vec<Vec<(usize, usize)>> = vec![];
// create a min-heap since our polynomials are sorted smallest to largest
let mut h: BinaryHeap<Reverse<_>> = BinaryHeap::with_capacity(self.nterms());
cache.insert(
Key {
index: 0,
monomials: &monomials,
},
vec![(0, 0)],
);
h.push(Reverse(Key {
index: 0,
monomials: &monomials,
}));
// i=merged_index[j] signifies that self[i]*other[j] has been merged
let mut merged_index = vec![0; rhs.nterms()];
// in_heap[j] signifies that other[j] is in the heap
let mut in_heap = vec![false; rhs.nterms()];
in_heap[0] = true;
while !h.is_empty() {
let cur_mon = h.pop().unwrap();
let mut coefficient = self.ring.zero();
let mut q = cache.remove(&cur_mon.0).unwrap();
for (i, j) in q.drain(..) {
self.ring.add_mul_assign(
&mut coefficient,
&self.coefficients[i],
&rhs.coefficients[j],
);
merged_index[j] = i + 1;
if i + 1 < self.nterms() && (j == 0 || merged_index[j - 1] > i + 1) {
let m = unsafe {
let b = &mut *monomials.get();
let index = b.1.len();
b.1.extend(
self.exponents(i + 1)
.iter()
.zip(rhs.exponents(j))
.map(|(e1, e2)| *e1 + *e2),
);
Key {
index,
monomials: &monomials,
}
};
if let Some(e) = cache.get_mut(&m) {
e.push((i + 1, j));
} else {
h.push(Reverse(m)); // only add when new
if let Some(mut qq) = q_cache.pop() {
qq.push((i + 1, j));
cache.insert(m, qq);
} else {
cache.insert(m, vec![(i + 1, j)]);
}
}
} else {
in_heap[j] = false;
}
if j + 1 < rhs.nterms() && !in_heap[j + 1] {
let m = unsafe {
let b = &mut *monomials.get();
let index = b.1.len();
b.1.extend(
self.exponents(i)
.iter()
.zip(rhs.exponents(j + 1))
.map(|(e1, e2)| *e1 + *e2),
);
Key {
index,
monomials: &monomials,
}
};
if let Some(e) = cache.get_mut(&m) {
e.push((i, j + 1));
} else {
h.push(Reverse(m)); // only add when new
if let Some(mut qq) = q_cache.pop() {
qq.push((i, j + 1));
cache.insert(m, qq);
} else {
cache.insert(m, vec![(i, j + 1)]);
}
}
in_heap[j + 1] = true;
}
}
q_cache.push(q);
if !self.ring.is_zero(&coefficient) {
res.coefficients.push(coefficient);
unsafe {
let b = &*monomials.get();
res.exponents
.extend_from_slice(&b.1[cur_mon.0.index..cur_mon.0.index + b.0]);
}
}
}
res
}
/// Heap multiplication, but with the exponents packed into a `u64`.
/// Each exponent is limited to 65535 if there are four or fewer variables,
/// or 255 if there are 8 or fewer variables.
fn heap_mul_packed_exp(
&self,
other: &MultivariatePolynomial<F, E, LexOrder>,
pack_u8: bool,
) -> MultivariatePolynomial<F, E, LexOrder> {
let mut res = self.zero_with_capacity(self.nterms().max(other.nterms()));
let pack_a: Vec<_> = if pack_u8 {
self.exponents_iter().map(|c| E::pack(c)).collect()
} else {
self.exponents_iter().map(|c| E::pack_u16(c)).collect()
};
let pack_b: Vec<_> = if pack_u8 {
other.exponents_iter().map(|c| E::pack(c)).collect()
} else {
other.exponents_iter().map(|c| E::pack_u16(c)).collect()
};
let mut cache: BTreeMap<u64, Vec<(usize, usize)>> = BTreeMap::new();
let mut q_cache: Vec<Vec<(usize, usize)>> = vec![];
// create a min-heap since our polynomials are sorted smallest to largest
let mut h: BinaryHeap<Reverse<u64>> = BinaryHeap::with_capacity(self.nterms());
let monom: u64 = pack_a[0] + pack_b[0];
cache.insert(monom, vec![(0, 0)]);
h.push(Reverse(monom));
// i=merged_index[j] signifies that self[i]*other[j] has been merged
let mut merged_index = vec![0; other.nterms()];
// in_heap[j] signifies that other[j] is in the heap
let mut in_heap = vec![false; other.nterms()];
in_heap[0] = true;
while let Some(cur_mon) = h.pop() {
let mut coefficient = self.ring.zero();
let mut q = cache.remove(&cur_mon.0).unwrap();
for (i, j) in q.drain(..) {
self.ring.add_mul_assign(
&mut coefficient,
&self.coefficients[i],
&other.coefficients[j],
);
merged_index[j] = i + 1;
if i + 1 < self.nterms() && (j == 0 || merged_index[j - 1] > i + 1) {
let m = pack_a[i + 1] + pack_b[j];
if let Some(e) = cache.get_mut(&m) {
e.push((i + 1, j));
} else {
h.push(Reverse(m)); // only add when new
if let Some(mut qq) = q_cache.pop() {
qq.push((i + 1, j));
cache.insert(m, qq);
} else {
cache.insert(m, vec![(i + 1, j)]);
}
}
} else {
in_heap[j] = false;
}
if j + 1 < other.nterms() && !in_heap[j + 1] {
let m = pack_a[i] + pack_b[j + 1];
if let Some(e) = cache.get_mut(&m) {
e.push((i, j + 1));
} else {
h.push(Reverse(m)); // only add when new
if let Some(mut qq) = q_cache.pop() {
qq.push((i, j + 1));
cache.insert(m, qq);
} else {
cache.insert(m, vec![(i, j + 1)]);
}
}
in_heap[j + 1] = true;
}
}
q_cache.push(q);
if !self.ring.is_zero(&coefficient) {
res.coefficients.push(coefficient);
let len = res.exponents.len();
res.exponents.resize(len + self.nvars(), E::zero());
if pack_u8 {
E::unpack(cur_mon.0, &mut res.exponents[len..len + self.nvars()]);
} else {
E::unpack_u16(cur_mon.0, &mut res.exponents[len..len + self.nvars()]);
}
}
}
res
}
/// Compute `self^pow` using a heap-based algorithm of "Sparse Polynomial Powering Using Heaps"
/// by Michael Monagan and Roman Pearce.
///
/// The caller must assure that the ring's cardinality is large enough to contain the exponents
/// after Kronecker mapping.
pub fn heap_pow(&self, pow: usize) -> Self {
if self.is_constant() {
return self.constant(self.ring.pow(&self.lcoeff(), pow as u64));
}
if self.coefficients.len() == 1 {
return self.monomial(
self.ring.pow(&self.coefficients[0], pow as u64),
self.exponents
.iter()
.map(|x| *x * E::from_i32(pow as i32))
.collect(),
);
}
#[inline(always)]
fn to_uni_var<E: Exponent>(s: &[E], max_degs_rev: &[usize]) -> Integer {
let mut shift = 1;
let mut res = Integer::from(s.last().unwrap().to_i32());
for (ee, &x) in s.iter().rev().skip(1).zip(max_degs_rev) {
shift *= x as u32;
res += ee.to_i32() as u32 * shift;
}
res
}
#[inline(always)]
fn from_uni_var<E: Exponent>(mut p: Integer, max_degs_rev: &[usize], exp: &mut [E]) {
for (ee, &x) in exp.iter_mut().rev().zip(max_degs_rev) {
*ee = E::from_i32(((&p % x as u64).to_i64().unwrap() as u32) as i32);
p /= x as u32;
}
}
let degree_bounds = (0..self.nvars())
.map(|v| self.degree_bounds(v))
.collect::<Vec<_>>();
let max_degs_rev = degree_bounds
.iter()
.rev()
.map(|v| (v.1 - v.0).to_i32() as usize * pow + 1)
.collect::<Vec<_>>();
let mut exp = vec![E::zero(); self.nvars()];
let mut f_exp: Vec<_> = self
.exponents_iter()
.map(|c| {
for ((ee, x), d) in exp.iter_mut().zip(c.iter()).zip(°ree_bounds) {
*ee = *x - d.0;
}
to_uni_var(&exp, &max_degs_rev)
})
.collect();
f_exp.reverse(); // descending order
let mut g_coeff = vec![self.ring.pow(self.coefficients.last().unwrap(), pow as u64)];
let mut g_exp = vec![f_exp[0].clone() * pow as u64];
let mut cache: BTreeMap<Integer, Vec<(usize, usize)>> = BTreeMap::new();
let mut q_cache: Vec<Vec<(usize, usize)>> = vec![];
// create a min-heap since our polynomials are sorted smallest to largest
let mut h: BinaryHeap<Integer> = BinaryHeap::with_capacity(self.nterms());
let monom = f_exp[1].clone() + &g_exp[0];
cache.insert(monom.clone(), vec![(1, 0)]);
h.push(monom);
// i=merged_index[j] signifies that self[i]*g[j] has been merged
let mut merged_index = vec![0; self.nterms()];
// in_heap[j] signifies that g[j] is in the heap
let mut in_heap = vec![false; self.nterms()];
in_heap[0] = true;
while let Some(cur_mon) = h.pop() {
let mut coefficient = self.ring.zero();
let mut q = cache.remove(&cur_mon).unwrap();
for (i, j) in q.drain(..) {
self.ring.add_mul_assign(
&mut coefficient,
&g_coeff[j],
&self.ring.mul(
self.coefficient_back(i),
&self
.ring
.nth(g_exp[j].clone() - f_exp[i].clone() * pow as u64),
),
);
if j + 1 >= merged_index.len() {
merged_index.resize(j + 2, 0);
in_heap.resize(j + 2, false);
}
merged_index[j] = i + 1;
if i + 1 < self.nterms() && (j == 0 || merged_index[j - 1] > i + 1) {
let m = f_exp[i + 1].clone() + &g_exp[j];
if let Some(e) = cache.get_mut(&m) {
e.push((i + 1, j));
} else {
h.push(m.clone()); // only add when new
if let Some(mut qq) = q_cache.pop() {
qq.push((i + 1, j));
cache.insert(m, qq);
} else {
cache.insert(m, vec![(i + 1, j)]);
}
}
} else {
in_heap[j] = false;
}
if j + 1 < g_exp.len() && !in_heap[j + 1] {
let m = f_exp[i].clone() + &g_exp[j + 1];
if let Some(e) = cache.get_mut(&m) {
e.push((i, j + 1));
} else {
h.push(m.clone()); // only add when new
if let Some(mut qq) = q_cache.pop() {
qq.push((i, j + 1));
cache.insert(m, qq);
} else {
cache.insert(m, vec![(i, j + 1)]);
}
}
in_heap[j + 1] = true;
}
}
q_cache.push(q);
if !self.ring.is_zero(&coefficient) {
g_exp.push(&cur_mon - &f_exp[0]);
let q = self
.ring
.try_div(
&coefficient,
&self.ring.mul(
self.coefficient_back(0),
&self.ring.nth(g_exp[0].clone() + &f_exp[0] - cur_mon),
),
)
.unwrap();
g_coeff.push(q);
if g_exp.len() >= in_heap.len() {
merged_index.resize(g_exp.len(), 0);
in_heap.resize(g_exp.len(), false);
}
if !in_heap[g_exp.len() - 1] {
let m = f_exp[1].clone() + &g_exp[g_exp.len() - 1];
if let Some(e) = cache.get_mut(&m) {
e.push((1, g_exp.len() - 1));
} else {
h.push(m.clone());
if let Some(mut qq) = q_cache.pop() {
qq.push((1, g_exp.len() - 1));
cache.insert(m, qq);
} else {
cache.insert(m, vec![(1, g_exp.len() - 1)]);
}
}
in_heap[g_exp.len() - 1] = true;
}
}
}
let mut res = self.zero();
for (c, e) in g_coeff.into_iter().zip(g_exp).rev() {
from_uni_var(e, &max_degs_rev, &mut exp);
for (ee, d) in exp.iter_mut().zip(°ree_bounds) {
*ee += d.0 * E::from_i32(pow as i32);
}
res.append_monomial(c, &exp);
}
res
}
}
impl<F: EuclideanDomain, E: PositiveExponent> MultivariatePolynomial<F, E, LexOrder> {
/// Convert the polynomial to one in a number field, where the variable
/// of the number field is moved into the coefficient.
pub fn to_number_field(
&self,
field: &AlgebraicExtension<F>,
) -> MultivariatePolynomial<AlgebraicExtension<F>, E> {
let var = &field.poly().get_vars_ref()[0];
let Some(var_index) = self.get_vars_ref().iter().position(|x| x == var) else {
return self.map_coeff(
|c| field.to_element(field.poly().constant(c.clone())),
field.clone(),
);
};
let polys = self.to_multivariate_polynomial_list(&[var_index], false);
// TODO: remove the variable from the variable map?
let mut poly =
MultivariatePolynomial::new(field, self.nterms().into(), self.variables.clone());
for (e, c) in polys {
let mut c2 = MultivariatePolynomial::new(
&self.ring,
c.nterms().into(),
Arc::new(vec![self.variables.as_ref()[var_index].clone()]),
);
c2.exponents = c
.exponents_iter()
.map(|x| x[var_index].to_i32() as u16)
.collect();
c2.coefficients = c.coefficients;
poly.append_monomial(field.to_element(c2), &e);
}
poly
}
/// Get the content from the coefficients.
pub fn content(&self) -> F::Element {
if self.coefficients.is_empty() {
return self.ring.zero();
}
let mut c = self.coefficients.first().unwrap().clone();
for cc in self.coefficients.iter().skip(1) {
// early return if possible (not possible for rationals)
if F::one_is_gcd_unit() && self.ring.is_one(&c) {
break;
}
c = self.ring.gcd(&c, cc);
}
c
}
/// Divide every coefficient with `other`.
pub fn div_coeff(mut self, other: &F::Element) -> Self {
for c in &mut self.coefficients {
let (quot, rem) = self.ring.quot_rem(c, other);
debug_assert!(self.ring.is_zero(&rem));
*c = quot;
}
self
}
/// Make the polynomial primitive by removing the content.
pub fn make_primitive(self) -> Self {
let c = self.content();
self.div_coeff(&c)
}
}
impl<F: EuclideanDomain, E: PositiveExponent> MultivariatePolynomial<F, E, LexOrder> {
/// Compute the remainder `self % div`.
pub fn rem(&self, div: &MultivariatePolynomial<F, E, LexOrder>) -> Self {
self.quot_rem(div, false).1
}
/// Divide two multivariate polynomials and return the quotient and remainder.
pub fn quot_rem(
&self,
div: &MultivariatePolynomial<F, E, LexOrder>,
abort_on_remainder: bool,
) -> (
MultivariatePolynomial<F, E, LexOrder>,
MultivariatePolynomial<F, E, LexOrder>,
) {
self.quot_rem_impl(div, abort_on_remainder)
}
/// Compute the p-adic expansion of the polynomial.
/// It returns `[a0, a1, a2, ...]` such that `a0 + a1 * p^1 + a2 * p^2 + ... = self`.
pub fn p_adic_expansion(&self, p: &Self) -> Vec<Self> {
if self.variables != p.variables {
let mut c1 = self.clone();
let mut c2 = p.clone();
c1.unify_variables(&mut c2);
return c1.p_adic_expansion(&c2);
}
let mut res = vec![];
let mut r = self.clone();
while !r.is_zero() {
let (q, rem) = r.quot_rem(p, false);
res.push(rem);
r = q;
}
res
}
}
impl<F: Ring, E: Exponent> MultivariatePolynomial<F, E, LexOrder> {
/// Divide `self` by `div` if there is no remainder, else return `None`.
pub fn try_div(
&self,
div: &MultivariatePolynomial<F, E, LexOrder>,
) -> Option<MultivariatePolynomial<F, E, LexOrder>> {
if div.is_zero() {
return None;
}
if self.variables != div.variables {
let mut c1 = self.clone();
let mut c2 = div.clone();
c1.unify_variables(&mut c2);
return c1.try_div(&c2);
}
if self.is_zero() {
return Some(self.clone());
}
// check if the leading coefficients divide
self.ring.try_div(&self.lcoeff(), &div.lcoeff())?;
if !self.is_polynomial() || !div.is_polynomial() {
// remove all negative exponents
let mut c1 = self.clone();
let mut c2 = div.clone();
let degrees = (0..self.nvars())
.map(|v| E::zero() - div.degree_bounds(v).0.min(E::zero()))
.collect::<Vec<_>>();
c1 = c1.mul_exp(°rees);
c2 = c2.mul_exp(°rees);
let mut degrees = (0..self.nvars())
.map(|v| E::zero() - self.degree_bounds(v).0.min(E::zero()))
.collect::<Vec<_>>();
c1 = c1.mul_exp(°rees);
let r = c1.try_div(&c2)?;
for d in &mut degrees {
*d = E::zero() - *d;
}
return Some(r.mul_exp(°rees));
}
if !self.is_polynomial() {
return None;
}
if (0..self.nvars()).any(|v| self.degree(v) < div.degree(v)) {
return None;
}
if self.ring.characteristic().is_zero() {
// test division of constant term (evaluation at x_i = 0)
let c = div.get_constant();
if !self.ring.is_zero(&c)
&& !self.ring.is_one(&c)
&& self.ring.try_div(&self.get_constant(), &c).is_none()
{
return None;
}
// test division at x_i = 1
let mut num = self.ring.zero();
for c in &self.coefficients {
self.ring.add_assign(&mut num, c);
}
let mut den = self.ring.zero();
for c in &div.coefficients {
self.ring.add_assign(&mut den, c);
}
if !self.ring.is_zero(&den)
&& !self.ring.is_one(&den)
&& self.ring.try_div(&num, &den).is_none()
{
return None;
}
}
let (a, b) = self.quot_rem_impl(div, true);
if b.nterms() == 0 { Some(a) } else { None }
}
/// Divide two multivariate polynomials and return the quotient and remainder.
///
/// The input must not have negative exponents.
fn quot_rem_impl(
&self,
div: &MultivariatePolynomial<F, E, LexOrder>,
abort_on_remainder: bool,
) -> (
MultivariatePolynomial<F, E, LexOrder>,
MultivariatePolynomial<F, E, LexOrder>,
) {
if div.is_zero() {
panic!("Cannot divide by 0 polynomial");
}
if self.is_zero() {
return (self.clone(), self.clone());
}
if div.is_one() {
return (self.clone(), self.zero());
}
if self.variables != div.variables {
let mut c1 = self.clone();
let mut c2 = div.clone();
c1.unify_variables(&mut c2);
return c1.quot_rem_impl(&c2, abort_on_remainder);
}
if self.nterms() == div.nterms() {
if self == div {
return (self.one(), self.zero());
}
// check if one is a multiple of the other
if let Some(q) = self.ring.try_div(&self.lcoeff(), &div.lcoeff())
&& self
.into_iter()
.zip(div)
.all(|(t1, t2)| t1.exponents == t2.exponents)
&& self
.into_iter()
.zip(div)
.all(|(t1, t2)| &self.ring.mul(t2.coefficient, &q) == t1.coefficient)
{
return (self.constant(q), self.zero());
}
}
if div.is_constant() {
let mut q = self.clone();
let dive = div.to_monomial_view(0);
if let Some(i) = div.ring.try_inv(dive.coefficient) {
return (q.mul_coeff(i), self.zero());
}
for c in &mut q.coefficients {
if let Some(quot) = self.ring.try_div(c, dive.coefficient) {
*c = quot;
} else {
return (self.zero(), self.clone());
}
}
return (q, self.zero());
}
// check if the division is univariate with the same variable
let degree_sum: Vec<_> = (0..self.nvars())
.map(|i| self.degree(i).to_i32() as usize + div.degree(i).to_i32() as usize)
.collect();
if div.ring.is_one(&div.lcoeff()) && degree_sum.iter().filter(|x| **x > 0).count() == 1 {
return self.quot_rem_univariate_monic(div);
}
let mut pack_u8 = true;
if self.nvars() <= 8
&& (0..self.nvars()).all(|i| {
let deg = self.degree(i).to_i32() as u32;
if deg > 127 {
pack_u8 = false;
}
deg <= 127 || self.nvars() <= 4 && deg <= 32767
})
{
self.heap_division_packed_exp(div, abort_on_remainder, pack_u8)
} else {
self.heap_division(div, abort_on_remainder)
}
}
/// Heap division for multivariate polynomials, using a cache so that only unique
/// monomial exponents appear in the heap.
/// Reference: "Sparse polynomial division using a heap" by Monagan, Pearce (2011)
///
/// The input must not have negative exponents.
fn heap_division(
&self,
div: &MultivariatePolynomial<F, E, LexOrder>,
abort_on_remainder: bool,
) -> (
MultivariatePolynomial<F, E, LexOrder>,
MultivariatePolynomial<F, E, LexOrder>,
) {
let mut q = self.zero_with_capacity(self.nterms());
let mut r = self.zero();
let mut div_monomial_in_heap = vec![false; div.nterms()];
let mut merged_index_of_div_monomial_in_quotient = vec![0; div.nterms()];
let mut cache: BTreeMap<Vec<E>, Vec<(usize, usize, bool)>> = BTreeMap::new();
let mut h: BinaryHeap<Vec<E>> = BinaryHeap::with_capacity(self.nterms());
let mut q_cache: Vec<Vec<(usize, usize, bool)>> = vec![];
let mut m = vec![E::zero(); div.nvars()];
let mut m_cache = vec![E::zero(); div.nvars()];
let mut c;
let mut k = 0;
while !h.is_empty() || k < self.nterms() {
if k < self.nterms() && (h.is_empty() || self.exponents_back(k) >= h.peek().unwrap()) {
for (s, e) in m.iter_mut().zip(self.exponents_back(k)) {
*s = *e;
}
c = self.coefficient_back(k).clone();
k += 1;
} else {
for (s, e) in m.iter_mut().zip(h.peek().unwrap().as_slice()) {
*s = *e;
}
c = self.ring.zero();
}
if let Some(monomial) = h.peek()
&& &m == monomial
{
h.pop().unwrap();
let mut qs = cache.remove(&m).unwrap();
for (i, j, next_in_divisor) in qs.drain(..) {
self.ring
.sub_mul_assign(&mut c, &q.coefficients[i], div.coefficient_back(j));
if next_in_divisor && j + 1 < div.nterms() {
// quotient heap product
for ((m, e1), e2) in m_cache
.iter_mut()
.zip(q.exponents(i))
.zip(div.exponents_back(j + 1))
{
*m = *e1 + *e2;
}
// TODO: make macro
if let Some(e) = cache.get_mut(&m_cache) {
e.push((i, j + 1, true));
} else {
h.push(m_cache.clone()); // only add when new
if let Some(mut qq) = q_cache.pop() {
qq.push((i, j + 1, true));
cache.insert(m_cache.clone(), qq);
} else {
cache.insert(m_cache.clone(), vec![(i, j + 1, true)]);
}
}
} else if !next_in_divisor {
merged_index_of_div_monomial_in_quotient[j] = i + 1;
if i + 1 < q.nterms()
&& (j == 1 // the divisor starts with the sub-leading term in the heap
|| merged_index_of_div_monomial_in_quotient[j - 1] > i + 1)
{
for ((m, e1), e2) in m_cache
.iter_mut()
.zip(q.exponents(i + 1))
.zip(div.exponents_back(j))
{
*m = *e1 + *e2;
}
if let Some(e) = cache.get_mut(&m_cache) {
e.push((i + 1, j, false));
} else {
h.push(m_cache.clone()); // only add when new
if let Some(mut qq) = q_cache.pop() {
qq.push((i + 1, j, false));
cache.insert(m_cache.clone(), qq);
} else {
cache.insert(m_cache.clone(), vec![(i + 1, j, false)]);
}
}
} else {
div_monomial_in_heap[j] = false;
}
if j + 1 < div.nterms() && !div_monomial_in_heap[j + 1] {
for ((m, e1), e2) in m_cache
.iter_mut()
.zip(q.exponents(i))
.zip(div.exponents_back(j + 1))
{
*m = *e1 + *e2;
}
if let Some(e) = cache.get_mut(&m_cache) {
e.push((i, j + 1, false));
} else {
h.push(m_cache.clone()); // only add when new
if let Some(mut qq) = q_cache.pop() {
qq.push((i, j + 1, false));
cache.insert(m_cache.clone(), qq);
} else {
cache.insert(m_cache.clone(), vec![(i, j + 1, false)]);
}
}
div_monomial_in_heap[j + 1] = true;
}
}
}
q_cache.push(qs);
}
if self.ring.is_zero(&c) {
continue;
}
if div.last_exponents().iter().zip(&m).all(|(ge, me)| me >= ge) {
if let Some(quot) = self.ring.try_div(&c, &div.lcoeff()) {
q.coefficients.push(quot);
} else if abort_on_remainder {
r = self.one();
return (q, r);
} else {
return (self.zero(), self.clone());
}
q.exponents.extend(
div.last_exponents()
.iter()
.zip(&m)
.map(|(ge, me)| *me - *ge),
);
if div.nterms() == 1 {
continue;
}
for ((m, e1), e2) in m_cache
.iter_mut()
.zip(q.last_exponents())
.zip(div.exponents_back(1))
{
*m = *e1 + *e2;
}
if q.nterms() < div.nterms() {
// using quotient heap
if let Some(e) = cache.get_mut(&m_cache) {
e.push((q.nterms() - 1, 1, true));
} else {
h.push(m_cache.clone()); // only add when new
if let Some(mut qq) = q_cache.pop() {
qq.push((q.nterms() - 1, 1, true));
cache.insert(m_cache.clone(), qq);
} else {
cache.insert(m_cache.clone(), vec![(q.nterms() - 1, 1, true)]);
}
}
} else if q.nterms() >= div.nterms() {
// using divisor heap
if !div_monomial_in_heap[1] {
div_monomial_in_heap[1] = true;
if let Some(e) = cache.get_mut(&m_cache) {
e.push((q.nterms() - 1, 1, false));
} else {
h.push(m_cache.clone()); // only add when new
if let Some(mut qq) = q_cache.pop() {
qq.push((q.nterms() - 1, 1, false));
cache.insert(m_cache.clone(), qq);
} else {
cache.insert(m_cache.clone(), vec![(q.nterms() - 1, 1, false)]);
}
}
}
} else {
// switch to divisor heap
for index in &mut merged_index_of_div_monomial_in_quotient {
*index = q.nterms() - 1;
}
debug_assert!(div_monomial_in_heap.iter().any(|c| !c));
div_monomial_in_heap[1] = true;
if let Some(e) = cache.get_mut(&m_cache) {
e.push((q.nterms() - 1, 1, false));
} else {
h.push(m_cache.clone()); // only add when new
if let Some(mut qq) = q_cache.pop() {
qq.push((q.nterms() - 1, 1, false));
cache.insert(m_cache.clone(), qq);
} else {
cache.insert(m_cache.clone(), vec![(q.nterms() - 1, 1, false)]);
}
}
}
} else if abort_on_remainder {
r = self.one();
return (q, r);
} else {
r.coefficients.push(c);
r.exponents.extend(&m);
}
}
// q and r have the highest monomials first
q.reverse_monomials();
r.reverse_monomials();
#[cfg(debug_assertions)]
{
if !(&q * div + r.clone() - self.clone()).is_zero() {
panic!("Division failed: ({self})/({div}): q={q}, r={r}");
}
}
(q, r)
}
/// Heap division, but with the exponents packed into a `u64`.
/// Each exponent is limited to 32767 if there are 5 or fewer variables,
/// or 127 if there are 8 or fewer variables, such that the last bit per byte can
/// be used to check for subtraction overflow, serving as a division test.
///
/// The input must not have negative exponents.
fn heap_division_packed_exp(
&self,
div: &MultivariatePolynomial<F, E, LexOrder>,
abort_on_remainder: bool,
pack_u8: bool,
) -> (
MultivariatePolynomial<F, E, LexOrder>,
MultivariatePolynomial<F, E, LexOrder>,
) {
let mut q = self.zero_with_capacity(self.nterms());
let mut r = self.zero();
let pack_a: Vec<_> = if pack_u8 {
self.exponents_iter().map(|c| E::pack(c)).collect()
} else {
self.exponents_iter().map(|c| E::pack_u16(c)).collect()
};
let pack_div: Vec<_> = if pack_u8 {
div.exponents_iter().map(|c| E::pack(c)).collect()
} else {
div.exponents_iter().map(|c| E::pack_u16(c)).collect()
};
let mut div_monomial_in_heap = vec![false; div.nterms()];
let mut merged_index_of_div_monomial_in_quotient = vec![0; div.nterms()];
let mut cache: BTreeMap<u64, Vec<(usize, usize, bool)>> = BTreeMap::new();
#[inline(always)]
fn divides(a: u64, b: u64, pack_u8: bool) -> Option<u64> {
let d = a.overflowing_sub(b).0;
if pack_u8 && (d & 9259542123273814144u64 == 0)
|| !pack_u8 && (d & 9223512776490647552u64 == 0)
{
Some(d)
} else {
None
}
}
let mut h: BinaryHeap<u64> = BinaryHeap::with_capacity(self.nterms());
let mut q_cache: Vec<Vec<(usize, usize, bool)>> = Vec::with_capacity(self.nterms());
let mut m;
let mut m_cache;
let mut c;
let mut q_exp = Vec::with_capacity(self.nterms());
let mut k = 0;
while !h.is_empty() || k < self.nterms() {
if k < self.nterms()
&& (h.is_empty() || pack_a[self.nterms() - k - 1] >= *h.peek().unwrap())
{
m = pack_a[self.nterms() - k - 1];
c = self.coefficient_back(k).clone();
k += 1;
} else {
m = *h.peek().unwrap();
c = self.ring.zero();
}
if let Some(monomial) = h.peek()
&& &m == monomial
{
h.pop().unwrap();
let mut qs = cache.remove(&m).unwrap();
for (i, j, next_in_divisor) in qs.drain(..) {
// TODO: use fraction-free routines
self.ring
.sub_mul_assign(&mut c, &q.coefficients[i], div.coefficient_back(j));
if next_in_divisor && j + 1 < div.nterms() {
// quotient heap product
m_cache = q_exp[i] + pack_div[div.nterms() - (j + 1) - 1];
// TODO: make macro
if let Some(e) = cache.get_mut(&m_cache) {
e.push((i, j + 1, true));
} else {
h.push(m_cache); // only add when new
if let Some(mut qq) = q_cache.pop() {
qq.push((i, j + 1, true));
cache.insert(m_cache, qq);
} else {
cache.insert(m_cache, vec![(i, j + 1, true)]);
}
}
} else if !next_in_divisor {
merged_index_of_div_monomial_in_quotient[j] = i + 1;
if i + 1 < q.nterms()
&& (j == 1 // the divisor starts with the sub-leading term in the heap
|| merged_index_of_div_monomial_in_quotient[j - 1] > i + 1)
{
m_cache = q_exp[i + 1] + pack_div[div.nterms() - j - 1];
if let Some(e) = cache.get_mut(&m_cache) {
e.push((i + 1, j, false));
} else {
h.push(m_cache); // only add when new
if let Some(mut qq) = q_cache.pop() {
qq.push((i + 1, j, false));
cache.insert(m_cache, qq);
} else {
cache.insert(m_cache, vec![(i + 1, j, false)]);
}
}
} else {
div_monomial_in_heap[j] = false;
}
if j + 1 < div.nterms() && !div_monomial_in_heap[j + 1] {
m_cache = q_exp[i] + pack_div[div.nterms() - (j + 1) - 1];
if let Some(e) = cache.get_mut(&m_cache) {
e.push((i, j + 1, false));
} else {
h.push(m_cache); // only add when new
if let Some(mut qq) = q_cache.pop() {
qq.push((i, j + 1, false));
cache.insert(m_cache, qq);
} else {
cache.insert(m_cache, vec![(i, j + 1, false)]);
}
}
div_monomial_in_heap[j + 1] = true;
}
}
}
q_cache.push(qs);
}
if self.ring.is_zero(&c) {
continue;
}
let q_e = divides(m, pack_div[pack_div.len() - 1], pack_u8);
if let Some(q_e) = q_e {
if let Some(quot) = self.ring.try_div(&c, &div.lcoeff()) {
q.coefficients.push(quot);
} else if abort_on_remainder {
r = self.one();
return (q, r);
} else {
return (self.zero(), self.clone());
}
let len = q.exponents.len();
q.exponents.resize(len + self.nvars(), E::zero());
if pack_u8 {
E::unpack(q_e, &mut q.exponents[len..len + self.nvars()]);
} else {
E::unpack_u16(q_e, &mut q.exponents[len..len + self.nvars()]);
}
q_exp.push(q_e);
if div.nterms() == 1 {
continue;
}
m_cache = q_exp.last().unwrap() + pack_div[pack_div.len() - 2];
if q.nterms() < div.nterms() {
// using quotient heap
if let Some(e) = cache.get_mut(&m_cache) {
e.push((q.nterms() - 1, 1, true));
} else {
h.push(m_cache); // only add when new
if let Some(mut qq) = q_cache.pop() {
qq.push((q.nterms() - 1, 1, true));
cache.insert(m_cache, qq);
} else {
cache.insert(m_cache, vec![(q.nterms() - 1, 1, true)]);
}
}
} else if q.nterms() >= div.nterms() {
// using divisor heap
if !div_monomial_in_heap[1] {
div_monomial_in_heap[1] = true;
if let Some(e) = cache.get_mut(&m_cache) {
e.push((q.nterms() - 1, 1, false));
} else {
h.push(m_cache); // only add when new
if let Some(mut qq) = q_cache.pop() {
qq.push((q.nterms() - 1, 1, false));
cache.insert(m_cache, qq);
} else {
cache.insert(m_cache, vec![(q.nterms() - 1, 1, false)]);
}
}
}
} else {
// switch to divisor heap
for index in &mut merged_index_of_div_monomial_in_quotient {
*index = q.nterms() - 1;
}
debug_assert!(div_monomial_in_heap.iter().any(|c| !c));
div_monomial_in_heap[1] = true;
if let Some(e) = cache.get_mut(&m_cache) {
e.push((q.nterms() - 1, 1, false));
} else {
h.push(m_cache); // only add when new
if let Some(mut qq) = q_cache.pop() {
qq.push((q.nterms() - 1, 1, false));
cache.insert(m_cache, qq);
} else {
cache.insert(m_cache, vec![(q.nterms() - 1, 1, false)]);
}
}
}
} else if abort_on_remainder {
r = self.one();
return (q, r);
} else {
r.coefficients.push(c);
let len = r.exponents.len();
r.exponents.resize(len + self.nvars(), E::zero());
if pack_u8 {
E::unpack(m, &mut r.exponents[len..len + self.nvars()]);
} else {
E::unpack_u16(m, &mut r.exponents[len..len + self.nvars()]);
}
}
}
// q and r have the highest monomials first
q.reverse_monomials();
r.reverse_monomials();
#[cfg(debug_assertions)]
{
if !(&q * div + r.clone() - self.clone()).is_zero() {
panic!("Division failed: ({self})/({div}): q={q}, r={r}");
}
}
(q, r)
}
}
impl<F: Field, E: Exponent, O: MonomialOrder> MultivariatePolynomial<F, E, O> {
/// Make the polynomial monic, i.e., make the leading coefficient `1` by
/// multiplying all monomials with `1/lcoeff`.
pub fn make_monic(self) -> Self {
if self.lcoeff() != self.ring.one() {
let ci = self.ring.inv(&self.lcoeff());
self.mul_coeff(ci)
} else {
self
}
}
}
impl<F: Field, E: PositiveExponent, O: MonomialOrder> MultivariatePolynomial<F, E, O> {
/// Integrate the polynomial w.r.t the variable `var`,
/// producing the antiderivative with zero constant.
pub fn integrate(&self, var: usize) -> Self {
debug_assert!(var < self.nvars());
if self.is_zero() {
return self.zero();
}
let mut res = self.zero_with_capacity(self.nterms());
let mut exp = vec![E::zero(); self.nvars()];
for x in self {
exp.copy_from_slice(x.exponents);
let pow = exp[var].to_u32() as u64;
exp[var] += E::one();
res.append_monomial(
self.ring
.div(x.coefficient, &self.ring.nth(Integer::from(pow) + 1)),
&exp,
);
}
res
}
}
impl<F: Field, E: PositiveExponent> MultivariatePolynomial<F, E, LexOrder> {
/// Optimized division routine for univariate polynomials over a field, which
/// makes the divisor monic first.
pub fn quot_rem_univariate(
&self,
div: &mut MultivariatePolynomial<F, E, LexOrder>,
) -> (
MultivariatePolynomial<F, E, LexOrder>,
MultivariatePolynomial<F, E, LexOrder>,
) {
if self.is_zero() {
return (self.clone(), self.clone());
}
if div.nterms() == 1 {
// calculate inverse once
let inv = self.ring.inv(&div.coefficients[0]);
if div.is_constant() {
let mut q = self.clone();
for c in &mut q.coefficients {
self.ring.mul_assign(c, &inv);
}
return (q, self.zero());
}
let mut q = self.zero_with_capacity(self.nterms());
let mut r = self.zero();
let dive = div.exponents(0);
for m in self.into_iter() {
if m.exponents.iter().zip(dive).all(|(a, b)| a >= b) {
q.coefficients.push(self.ring.mul(m.coefficient, &inv));
for (ee, ed) in m.exponents.iter().zip(dive) {
q.exponents.push(*ee - *ed);
}
} else {
r.coefficients.push(m.coefficient.clone());
r.exponents.extend(m.exponents);
}
}
return (q, r);
}
// normalize the lcoeff to 1 to prevent a costly inversion
if !self.ring.is_one(&div.lcoeff()) {
let o = div.lcoeff();
let inv = self.ring.inv(&div.lcoeff());
for c in &mut div.coefficients {
self.ring.mul_assign(c, &inv);
}
let mut res = self.quot_rem_univariate_monic(div);
for c in &mut res.0.coefficients {
self.ring.mul_assign(c, &inv);
}
for c in &mut div.coefficients {
self.ring.mul_assign(c, &o);
}
return res;
}
self.quot_rem_univariate_monic(div)
}
/// Compute self^n % m where m is a polynomial
pub fn exp_mod_univariate(&self, mut n: Integer, m: &mut Self) -> Self {
if n.is_zero() {
return self.one();
}
// use binary exponentiation and mod at every stage
let mut x = self.quot_rem_univariate(m).1;
let mut y = self.one();
while !n.is_one() {
if (&n % &Integer::Single(2)).is_one() {
y = (&y * &x).quot_rem_univariate(m).1;
n -= &Integer::one();
}
x = (&x * &x).quot_rem_univariate(m).1;
n /= 2;
}
(x * &y).quot_rem_univariate(m).1
}
/// Perform a fast univariate division of `self` by `div`, using a precomputed inverse of `div` reversed.
pub fn quot_rem_univariate_fast(&self, div: &Self, var: usize, inv_div: &Self) -> (Self, Self) {
if div.is_zero() {
panic!("Cannot divide by 0 polynomial");
}
if self.is_zero() {
return (self.clone(), self.clone());
}
let deg_a = *self.last_exponents().iter().max().unwrap();
let deg_b = *div.last_exponents().iter().max().unwrap();
if deg_a < deg_b {
return (self.zero(), self.clone());
}
let m = deg_a - deg_b;
let mut self_i = self.clone();
self_i.reverse();
let mut q = (&self_i * &inv_div.mod_var(var, m + E::one())).mod_var(var, m + E::one());
q.reverse();
if q.degree(var) < m {
let mut exp = vec![E::zero(); self.nvars()];
exp[var] = m - q.degree(var);
q = q.mul_exp(&exp);
}
let r = self - &(div * &q);
(q, r)
}
/// Compute `self^n % m` where `m` is a polynomial, using a precomputed inverse of `m` reversed.
pub fn exp_mod_univariate_fast(
&mut self,
var: usize,
mut n: Integer,
m: &Self,
m_inv: &Self,
) -> Self {
if n.is_zero() {
return self.one();
}
// use binary exponentiation and mod at every stage
let mut x = self.quot_rem_univariate_fast(m, var, m_inv).1;
let mut y = self.one();
while !n.is_one() {
if (&n % &Integer::Single(2)).is_one() {
y = (&y * &x).quot_rem_univariate_fast(m, var, m_inv).1;
n -= &Integer::one();
}
x = (&x * &x).quot_rem_univariate_fast(m, var, m_inv).1;
n /= 2;
}
(x * &y).quot_rem_univariate_fast(m, var, m_inv).1
}
/// Compute `(g, s, t)` where `self * s + other * t = g`
/// by means of the extended Euclidean algorithm.
pub fn eea_univariate(&self, other: &Self) -> (Self, Self, Self) {
let mut r0 = self.clone().make_monic();
let mut r1 = other.clone().make_monic();
let mut s0 = self.constant(self.ring.inv(&self.lcoeff()));
let mut s1 = self.zero();
let mut t0 = self.zero();
let mut t1 = self.constant(self.ring.inv(&other.lcoeff()));
while !r1.is_zero() {
let (q, r) = r0.quot_rem_univariate(&mut r1);
if self.ring.is_zero(&r.lcoeff()) {
return (r1, s1, t1);
}
let a = self.ring.inv(&r.lcoeff());
(r1, r0) = (r.mul_coeff(a.clone()), r1);
(s1, s0) = ((s0 - &q * &s1).mul_coeff(a.clone()), s1);
(t1, t0) = ((t0 - q * &t1).mul_coeff(a), t1);
}
(r0, s0, t0)
}
/// Compute `(s1,...,n2)` where `A0 * s0 + ... + An * sn = g`
/// where `Ai = prod(polys[j], j != i)`
/// by means of the extended Euclidean algorithm.
///
/// The `polys` must be pairwise co-prime.
pub fn diophantine_univariate(polys: &mut [Self], b: &Self) -> Vec<Self> {
let mut cur = polys.last().unwrap().clone();
let mut a = vec![cur.clone()];
for x in polys[1..].iter().rev().skip(1) {
cur = cur * x;
a.push(cur.clone());
}
a.reverse();
let mut ss = vec![];
let mut cur_s = b.clone();
for (p, aa) in polys.iter_mut().zip(&mut a) {
let (g, s, t) = p.eea_univariate(aa);
debug_assert!(g.is_one());
let new_s = (t * &cur_s).quot_rem_univariate(p).1;
ss.push(new_s);
cur_s = (s * &cur_s).quot_rem_univariate(aa).1;
}
ss.push(cur_s);
ss
}
/// Find a rational fraction `n(x)/d(x)`, the Pade approximant,
/// such that `d(x)*self-n(x)=0 mod x^(deg_n+deg_d+1)` and
/// `deg(d(x)) <= deg_d` and `deg(n(x) <= deg_n` using the extended Euclidean algorithm.
pub fn rational_approximant_univariate(&self, deg_n: u32, deg_d: u32) -> Option<(Self, Self)>
where
F: PolynomialGCD<E>,
{
let Some(var) = self.last_exponents().iter().position(|x| *x > E::zero()) else {
return Some((self.clone(), self.one()));
};
let mut exp = self.last_exponents().to_vec();
exp[var] = E::from_u32(deg_n) + E::from_u32(deg_d) + E::one();
let mut v0 = self.monomial(self.ring.one(), exp);
let mut v1 = self.zero();
let mut w0 = self.clone();
let mut w1 = self.one();
while w0.degree(var).to_u32() > deg_n {
let (q, r) = v0.quot_rem_univariate(&mut w0);
(w1, v1) = (v1 - q * &w1, w1);
(v0, w0) = (w0, r);
}
// TODO: normalize denominator?
let r = w0.gcd(&w1);
Some((w0 / &r, w1 / &r))
}
/// Shift a variable `var` to `var+shift`, using an optimized routine that
/// uses a power cache. If working in a finite field, the characteristic
/// should be larger than the degree of the polynomial.
pub fn shift_var_cached(&self, var: usize, shift: &F::Element) -> Self {
let d = self.degree(var).to_u32() as usize;
let y_poly = self.to_univariate_polynomial_list(var);
let mut sample_powers = Vec::with_capacity(d + 1);
let mut accum = self.ring.one();
sample_powers.push(self.ring.one());
for _ in 0..d {
self.ring.mul_assign(&mut accum, shift);
sample_powers.push(accum.clone());
}
let mut v = vec![self.zero(); d + 1];
for (x_poly, p) in y_poly {
let i = p.to_u32() as usize;
v[i] = x_poly.mul_coeff(sample_powers[i].clone());
}
for k in 0..d {
for j in (k..d).rev() {
v[j] = &v[j] + &v[j + 1];
}
}
let mut accum_inv = self.ring.one();
let sample_point_inv = self.ring.inv(shift);
let mut shifted_coefficients = Vec::with_capacity(v.len());
for (i, mut v) in v.into_iter().enumerate() {
v = v.mul_coeff(accum_inv.clone());
for x in v.exponents.chunks_mut(self.nvars()) {
x[var] = E::from_u32(i as u32);
}
shifted_coefficients.push(v);
self.ring.mul_assign(&mut accum_inv, &sample_point_inv);
}
self.merge_shifted_univariate_coefficients(shifted_coefficients)
}
}
impl<R: Ring, E: Exponent> Derivable for PolynomialRing<R, E> {
type Variable = PolyVariable;
fn derivative(
&self,
p: &MultivariatePolynomial<R, E>,
x: &PolyVariable,
) -> MultivariatePolynomial<R, E> {
if let Some(pos) = p.get_vars_ref().iter().position(|v| v == x) {
p.derivative(pos)
} else {
self.zero()
}
}
}
impl<R: EuclideanDomain, E: Exponent> MultivariatePolynomial<AlgebraicExtension<R>, E> {
/// Convert the polynomial to a multivariate polynomial that contains the
/// variable in the number field.
pub fn from_number_field(&self) -> MultivariatePolynomial<R, E> {
let var = &self.ring.poly().get_vars_ref()[0];
let (var_map, var_index) = if let Some(p) =
self.get_vars_ref().iter().position(|v| v == var)
{
if self.degree(p) > E::zero() {
panic!(
"The variable of the minimal polynomial of the coefficient field also appears in the polynomial"
);
}
(self.variables.clone(), p)
} else {
let p = self.get_vars_ref().len();
let mut v = self.get_vars_ref().to_vec();
v.push(var.clone());
(Arc::new(v), p)
};
let mut poly =
MultivariatePolynomial::new(&self.ring.poly().ring, self.nterms().into(), var_map);
let mut exp = vec![E::zero(); poly.nvars()];
for t in self {
exp[..self.nvars()].copy_from_slice(t.exponents);
for t2 in &t.coefficient.poly {
exp[var_index] = E::from_i32(t2.exponents[0].to_i32());
poly.append_monomial(t2.coefficient.clone(), &exp);
}
}
poly
}
}
impl<E: Exponent> MultivariatePolynomial<IntegerRing, E> {
/// Compute the polynomial that is congruent to `self` modulo `m` and `other` modulo `p` using the Chinese Remainder Theorem.
pub fn chinese_remainder(&self, other: &Self, m: &Integer, p: &Integer) -> Self {
let mut i = 0;
let mut j = 0;
let mut res = self.zero();
while i < self.nterms() || j < other.nterms() {
let (exp, mut c1, mut c2) = if i < self.nterms() && j < other.nterms() {
match self.exponents(i).cmp(other.exponents(j)) {
std::cmp::Ordering::Equal => {
i += 1;
j += 1;
(
self.exponents(i - 1),
self.coefficients[i - 1].clone(),
other.coefficients[j - 1].clone(),
)
}
std::cmp::Ordering::Less => {
i += 1;
(
self.exponents(i - 1),
self.coefficients[i - 1].clone(),
0.into(),
)
}
std::cmp::Ordering::Greater => {
j += 1;
(
other.exponents(j - 1),
0.into(),
other.coefficients[j - 1].clone(),
)
}
}
} else if i < self.nterms() {
i += 1;
(
self.exponents(i - 1),
self.coefficients[i - 1].clone(),
0.into(),
)
} else {
j += 1;
(
other.exponents(j - 1),
0.into(),
other.coefficients[j - 1].clone(),
)
};
if c1.is_negative() {
c1 += m;
}
if c2.is_negative() {
c2 += p;
}
let coeff = Integer::chinese_remainder(c1, c2, m.clone(), p.clone());
res.append_monomial(coeff, exp);
}
res
}
}
impl<E: Exponent> From<&MultivariatePolynomial<IntegerRing, E>>
for MultivariatePolynomial<RationalField, E>
{
fn from(val: &MultivariatePolynomial<IntegerRing, E>) -> Self {
MultivariatePolynomial {
coefficients: val.coefficients.iter().map(|x| x.into()).collect(),
exponents: val.exponents.clone(),
ring: Q,
variables: val.variables.clone(),
_phantom: PhantomData,
}
}
}
/// View object for a term in a multivariate polynomial.
#[derive(Copy, Clone, Debug)]
pub struct MonomialView<'a, F: 'a + Ring, E: 'a + Exponent> {
pub coefficient: &'a F::Element,
pub exponents: &'a [E],
}
/// Iterator over terms in a multivariate polynomial.
pub struct MonomialViewIterator<'a, F: Ring, E: Exponent, O: MonomialOrder> {
poly: &'a MultivariatePolynomial<F, E, O>,
index: usize,
}
impl<'a, F: Ring, E: Exponent, O: MonomialOrder> Iterator for MonomialViewIterator<'a, F, E, O> {
type Item = MonomialView<'a, F, E>;
#[inline]
fn next(&mut self) -> Option<Self::Item> {
if self.index == self.poly.nterms() {
None
} else {
let view = MonomialView {
coefficient: &self.poly.coefficients[self.index],
exponents: self.poly.exponents(self.index),
};
self.index += 1;
Some(view)
}
}
}
impl<'a, F: Ring, E: Exponent, O: MonomialOrder> IntoIterator
for &'a MultivariatePolynomial<F, E, O>
{
type Item = MonomialView<'a, F, E>;
type IntoIter = MonomialViewIterator<'a, F, E, O>;
#[inline]
fn into_iter(self) -> Self::IntoIter {
Self::IntoIter {
poly: self,
index: 0,
}
}
}
// General implementation for F::Element blocked on https://github.com/rust-lang/rust/issues/20400
impl<E: Exponent, O: MonomialOrder> PartialEq<Integer>
for MultivariatePolynomial<IntegerRing, E, O>
{
#[inline]
fn eq(&self, other: &Integer) -> bool {
self.is_constant() && self.get_constant() == *other
}
}
impl<R: Ring + FractionNormalization + EuclideanDomain, E: Exponent, O: MonomialOrder>
PartialEq<Fraction<R>> for MultivariatePolynomial<FractionField<R>, E, O>
{
#[inline]
fn eq(&self, other: &Fraction<R>) -> bool {
self.is_constant() && self.get_constant() == *other
}
}
#[cfg(test)]
mod test {
use crate::{atom::AtomCore, domains::integer::Z, domains::rational::Q, parse, symbol};
#[test]
fn mul_packed() {
let p1 = parse!("v1^2+v2^3*v3*+3*v1^4+4*v2*v3").to_polynomial::<_, u8>(&Z, None);
let b = &p1 * &p1;
let r = parse!(
"16*v2^2*v3^2+8*v1^2*v2*v3+v1^4+24*v1^4*v2^4*v3^2+6*v1^6*v2^3*v3+9*v1^8*v2^6*v3^2"
);
assert_eq!(b.to_expression(), r)
}
#[test]
fn mul_full() {
let p1 = parse!("v1^2+v2^3*v3*+3*v1^4+4*v2*v3+v4+v5+v6*v1*v2+v7*v5+v8+v9*v8")
.to_polynomial::<_, u8>(&Z, None);
let b = &p1 * &p1;
let r = parse!(
"16*v2^2*v3^2+8*v1*v2^2*v3*v6+8*v1^2*v2*v3+v1^2*v2^2*v6^2+2*v1^3*v2*v6+v1^4+24*v1^4*v2^4*v3^2+6*v1^5*v2^4*v3*v6+6*v1^6*v2^3*v3+9*v1^8*v2^6*v3^2+8*v8*v2*v3+8*v8*v2*v3*v9+2*v8*v1*v2*v6+2*v8*v1*v2*v9*v6+2*v8*v1^2+2*v8*v1^2*v9+6*v8*v1^4*v2^3*v3+6*v8*v1^4*v2^3*v3*v9+v8^2+2*v8^2*v9+v8^2*v9^2+8*v5*v2*v3+8*v5*v2*v3*v7+2*v5*v1*v2*v6+2*v5*v1*v2*v7*v6+2*v5*v1^2+2*v5*v1^2*v7+6*v5*v1^4*v2^3*v3+6*v5*v1^4*v2^3*v3*v7+2*v5*v8+2*v5*v8*v9+2*v5*v8*v7+2*v5*v8*v7*v9+v5^2+2*v5^2*v7+v5^2*v7^2+8*v4*v2*v3+2*v4*v1*v2*v6+2*v4*v1^2+6*v4*v1^4*v2^3*v3+2*v4*v8+2*v4*v8*v9+2*v4*v5+2*v4*v5*v7+v4^2"
);
assert_eq!(b.to_expression(), r)
}
#[test]
fn div_packed() {
let p1 = parse!("(v1+v2*5+v3*v2+v1*v2*v3)(v1+v2+v3)").to_polynomial::<_, u8>(&Z, None);
let p2 = parse!("v1+v2+v3+1").to_polynomial::<_, u8>(&Z, p1.variables.clone());
let (q, r) = p1.quot_rem(&p2, false);
assert_eq!(q.to_expression(), parse!("-1+5*v2+v1+v1*v2*v3"));
assert_eq!(r.to_expression(), parse!("1+v3-4*v2+v2*v3^2+v2^2*v3"));
}
#[test]
fn div_full() {
let p1 = parse!("(v1+v2*5+v3*v2+v1*v2*v3+v4+v5+v6+v7+v8+v9*v8)(v1+v2+v3)")
.to_polynomial::<_, u8>(&Z, None);
let p2 = parse!("v1+v2+v3+1").to_polynomial::<_, u8>(&Z, p1.variables.clone());
let (q, r) = p1.quot_rem(&p2, false);
assert_eq!(
q.to_expression(),
parse!("-1+v8+v8*v9+v7+v6+v5+v4+5*v2+v1+v1*v2*v3")
);
assert_eq!(
r.to_expression(),
parse!("1-v8-v8*v9-v7-v6-v5-v4+v3-4*v2+v2*v3^2+v2^2*v3")
);
}
#[test]
fn fuse_variables() {
let p1 = parse!("v1+v2").to_polynomial::<_, u8>(&Z, None);
let p2 = parse!("v4").to_polynomial::<_, u8>(&Z, None);
let p3 = parse!("v3").to_polynomial::<_, u8>(&Z, p1.variables.clone());
let r = p1 * &p2 + p3;
assert_eq!(
r.get_vars_ref(),
&[symbol!("v1"), symbol!("v2"), symbol!("v4"), symbol!("v3")]
);
}
#[test]
fn fast_exp_mod() {
let mut a = parse!("x^3 + 2*x + 3").to_polynomial::<_, u8>(&Q, None);
let b = parse!("x^2 + x + 2").to_polynomial::<_, u8>(&Q, None);
let mut b_rev = b.clone();
b_rev.reverse();
let c = b_rev.inverse_univariate(0, 2.into());
let d = a.exp_mod_univariate_fast(0, 4.into(), &b, &c);
assert_eq!(d.to_expression(), parse!("367+333*x"));
}
}