//! Factorization methods for multivariate polynomials
//! that implement [Factorize].
use std::{borrow::Cow, cmp::Reverse};
use ahash::{HashMap, HashSet, HashSetExt};
use rand::{Rng, SeedableRng, rng, rngs::StdRng};
use tracing::debug;
use crate::{
GLOBAL_SETTINGS,
combinatorics::CombinationIterator,
domains::{
EuclideanDomain, Field, InternalOrdering, Ring, RingOps, Set,
algebraic_number::{AlgebraicExtension, GaloisField},
finite_field::{
FiniteField, FiniteFieldCore, FiniteFieldWorkspace, PrimeIteratorU64, ToFiniteField,
Zp, Zp64,
},
integer::{Integer, IntegerRing, Z, gcd_unsigned},
rational::{Q, RationalField},
},
};
use super::{
LexOrder, PositiveExponent,
gcd::{MAX_RNG_PREFACTOR, POW_CACHE_SIZE, PolynomialGCD},
polynomial::MultivariatePolynomial,
};
const SPARSE_MDP_SAMPLE_BASE_ATTEMPTS: usize = 1;
// Keep the old bivariate start for cases where specializing to two variables
// does not collapse too many terms onto the same bivariate support.
const INTEGER_FACTOR_BIVARIATE_SPARSE_BOX_DENSITY_THRESHOLD: f64 = 3.0;
#[derive(Clone, Copy, Debug, Eq, PartialEq)]
enum IntegerFactorStart {
Auto,
Univariate,
Bivariate,
Disabled,
}
fn integer_factor_start_mode() -> IntegerFactorStart {
match (
GLOBAL_SETTINGS
.use_univariate_factorization
.load(std::sync::atomic::Ordering::Relaxed),
GLOBAL_SETTINGS
.use_bivariate_factorization
.load(std::sync::atomic::Ordering::Relaxed),
) {
(true, true) => IntegerFactorStart::Auto,
(true, false) => IntegerFactorStart::Univariate,
(false, true) => IntegerFactorStart::Bivariate,
(false, false) => IntegerFactorStart::Disabled,
}
}
/// A polynomial that can be factorized.
pub trait Factorize: Sized {
/// Perform a square-free factorization.
/// The output is `a_1^e1*...*a_n^e_n`
/// where each `a_i` is relative prime.
fn square_free_factorization(&self) -> Vec<(Self, usize)>;
/// Factor a polynomial over its coefficient ring.
fn factor(&self) -> Vec<(Self, usize)>;
fn is_irreducible(&self) -> bool;
}
impl<F: EuclideanDomain + PolynomialGCD<E>, E: PositiveExponent>
MultivariatePolynomial<F, E, LexOrder>
{
/// Find factors that do not contain all variables.
pub fn factor_separable(&self) -> Vec<Self> {
let mut stripped = self.clone();
let mut factors = vec![];
for x in 0..self.nvars() {
if stripped.degree(x) == E::zero() {
continue;
}
let c = stripped.to_univariate_polynomial_list(x);
let cs = c.into_iter().map(|x| x.0).collect();
let gcd = PolynomialGCD::gcd_multiple(cs);
if !gcd.is_constant() {
stripped = stripped / &gcd;
let mut fs = gcd.factor_separable();
factors.append(&mut fs);
}
}
factors.push(stripped);
factors
}
/// Perform a square free factorization using Yun's algorithm.
///
/// The characteristic of the ring must be 0 and all variables
/// must occur in every factor.
fn square_free_factorization_0_char(&self) -> Vec<(Self, usize)> {
if self.is_constant() {
if self.is_one() {
return vec![];
} else {
return vec![(self.clone(), 1)];
}
}
// any variable can be selected
// select the one with the lowest degree
let lowest_rank_var = (0..self.nvars())
.filter_map(|x| {
let d = self.degree(x);
if d > E::zero() { Some((x, d)) } else { None }
})
.min_by_key(|a| a.1)
.unwrap()
.0;
let b = self.derivative(lowest_rank_var);
let c = self.gcd(&b);
if c.is_one() {
return vec![(self.clone(), 1)];
}
let mut factors = vec![];
let mut w = self / &c;
let mut y = &b / &c;
let mut i = 1;
while !w.is_constant() {
let z = y - w.derivative(lowest_rank_var);
let g = w.gcd(&z);
w = w / &g;
y = z / &g;
if !g.is_one() {
factors.push((g, i));
}
i += 1
}
factors
}
/// Use Newton's polygon method to test if a bivariate polynomial is irreducible.
/// If this method returns `false`, the test is inconclusive.
///
/// The polynomial must have overall factors of single variables removed.
fn bivariate_irreducibility_test(&self) -> bool {
/// Compute the convex hull via the Monotone chain algorithm.
fn convex_hull(mut points: Vec<(isize, isize)>) -> Vec<(isize, isize)> {
points.sort();
if points.len() < 2 {
return points;
}
// Cross product of o-a and o-b vectors, positive means ccw turn, negative means cw turn and 0 means collinear.
fn cross(o: &(isize, isize), a: &(isize, isize), b: &(isize, isize)) -> isize {
(a.0 - o.0) * (b.1 - o.1) - (a.1 - o.1) * (b.0 - o.0)
}
let mut lower = vec![];
let mut upper = vec![];
for (t, rev) in [(&mut lower, false), (&mut upper, true)] {
for i in 0..points.len() {
let p = if rev {
points[points.len() - 1 - i]
} else {
points[i]
};
while t.len() >= 2 && cross(&t[t.len() - 2], &t[t.len() - 1], &p) <= 0 {
t.pop();
}
t.push(p);
}
}
lower.pop();
upper.pop();
lower.extend(upper);
lower
}
let vars: Vec<_> = (0..self.nvars())
.filter(|v| self.degree(*v) > E::zero())
.collect();
if vars.len() != 2 {
return false;
}
let points = self
.exponents
.chunks(self.nvars())
.map(|e| (e[vars[0]].to_u32() as isize, e[vars[1]].to_u32() as isize))
.collect();
let hull = convex_hull(points);
match hull.len() {
2 => {
let x_deg = hull[0].0.abs_diff(hull[1].0);
let y_deg = hull[0].1.abs_diff(hull[1].1);
gcd_unsigned(x_deg as u64, y_deg as u64) == 1
}
3 => {
// the hull has the form (n, 0), (0, m), (u, v)
let (mut n, mut m, mut u, mut v) = (-1, -1, -1, -1);
for (x, y) in hull {
if x != 0 && y == 0 {
n = x;
} else if y != 0 && x == 0 {
m = y;
} else {
u = x;
v = y;
}
}
n != -1
&& m != -1
&& u != -1
&& v != -1
&& gcd_unsigned(
gcd_unsigned(gcd_unsigned(n as u64, m as u64), u as u64),
v as u64,
) == 1
}
_ => false,
}
}
}
impl<R: EuclideanDomain, E: PositiveExponent> MultivariatePolynomial<R, E, LexOrder> {
/// Check if a parse lift is possible.
#[allow(dead_code)]
fn sparse_lift_possible(&self, factors: &[Self], order: &[usize]) -> bool {
// check if all bivariate monomials occur in the product of factors
let mut all_monomials = HashSet::with_capacity(self.nterms());
for e in self.exponents.chunks(self.nvars()) {
all_monomials.insert((e[order[0]], e[order[1]]));
}
let mut total = factors[0].clone();
for f in &factors[1..] {
total = &total * f;
}
let mut all_monomials_in_factors = HashSet::with_capacity(self.nterms());
for e in total.exponents.chunks(total.nvars()) {
all_monomials_in_factors.insert((e[order[0]], e[order[1]]));
}
all_monomials == all_monomials_in_factors
}
}
impl<E: PositiveExponent> Factorize for MultivariatePolynomial<IntegerRing, E, LexOrder> {
fn square_free_factorization(&self) -> Vec<(Self, usize)> {
if self.is_zero() {
return vec![];
}
let mut c = self.content();
let stripped = self.clone().div_coeff(&c);
let mut factors = vec![];
let fs = stripped.factor_separable();
for mut f in fs {
// make sure f is primitive
if f.lcoeff().is_negative() {
c = -c;
f = -f;
}
let mut nf = f.square_free_factorization_0_char();
factors.append(&mut nf);
}
if !c.is_one() {
factors.insert(0, (self.constant(c), 1));
}
if factors.is_empty() {
factors.push((self.one(), 1))
}
factors
}
fn factor(&self) -> Vec<(Self, usize)> {
let sf = self.square_free_factorization();
let mut factors = vec![];
let mut degrees = vec![0; self.nvars()];
for (f, p) in sf {
debug!("SFF {} {}", f, p);
let mut var_count = 0;
for (v, d) in degrees.iter_mut().enumerate() {
*d = f.degree(v).to_u32() as usize;
if *d > 0 {
var_count += 1;
}
}
match var_count {
0 | 1 => factors.extend(f.factor_reconstruct().into_iter().map(|ff| (ff, p))),
2 => {
let mut order: Vec<_> = degrees
.iter()
.enumerate()
.filter(|(_, d)| **d > 0)
.collect();
order.sort_by_key(|o| Reverse(o.1));
let order: Vec<_> = order.into_iter().map(|(v, _)| v).collect();
factors.extend(
f.bivariate_factor_reconstruct(order[0], order[1])
.into_iter()
.map(|ff| (ff, p)),
)
}
_ => {
// select the variable with the smallest leading coefficient and the highest degree to be first
let mut lcoeff_length = vec![0; self.nvars()];
for x in f.exponents_iter() {
for ((lc, e), d) in lcoeff_length.iter_mut().zip(x).zip(°rees) {
if e.to_i32() as usize == *d {
*lc += 1;
}
}
}
let first = (0..self.nvars())
.min_by(|a, b| {
lcoeff_length[*a]
.cmp(&lcoeff_length[*b])
.then_with(|| degrees[*b].cmp(°rees[*a]))
})
.unwrap();
// TODO: find better order
let mut order: Vec<_> = degrees
.iter()
.enumerate()
.filter(|(_, d)| **d > 0)
.collect();
order.sort_by_key(|o| {
if o.0 == first {
Reverse(&usize::MAX)
} else {
Reverse(o.1)
}
});
let mut order: Vec<_> = order.into_iter().map(|(v, _)| v).collect();
factors.extend(
f.multivariate_factorization(&mut order, 10, None)
.into_iter()
.map(|ff| (ff, p)),
)
}
}
}
factors
}
fn is_irreducible(&self) -> bool {
let mut sf = self.square_free_factorization();
if sf.len() > 1 {
return false;
}
let (f, _) = sf.pop().unwrap();
let mut degrees = vec![0; self.nvars()];
let mut var_count = 0;
for (v, d) in degrees.iter_mut().enumerate() {
*d = f.degree(v).to_u32() as usize;
if *d > 0 {
var_count += 1;
}
}
match var_count {
0 | 1 => f.factor_reconstruct().len() == 1,
2 => {
let mut order: Vec<_> = degrees
.iter()
.enumerate()
.filter(|(_, d)| **d > 0)
.collect();
order.sort_by_key(|o| Reverse(o.1));
let order: Vec<_> = order.into_iter().map(|(v, _)| v).collect();
f.bivariate_factor_reconstruct(order[0], order[1]).len() == 1
}
_ => {
// TODO: find better order
let mut order: Vec<_> = degrees
.iter()
.enumerate()
.filter(|(_, d)| **d > 0)
.collect();
order.sort_by_key(|o| Reverse(o.1));
let mut order: Vec<_> = order.into_iter().map(|(v, _)| v).collect();
f.multivariate_factorization(&mut order, 10, None).len() == 1
}
}
}
}
impl<E: PositiveExponent> Factorize for MultivariatePolynomial<RationalField, E, LexOrder> {
fn square_free_factorization(&self) -> Vec<(Self, usize)> {
if self.is_zero() {
return vec![];
}
let c = self.content();
let stripped = self.map_coeff(
|coeff| {
let coeff = self.ring.div(coeff, &c);
debug_assert!(coeff.is_integer());
coeff.numerator()
},
Z,
);
let fs = stripped.square_free_factorization();
let mut factors: Vec<_> = fs
.into_iter()
.map(|(f, e)| (f.map_coeff(|coeff| coeff.into(), Q), e))
.collect();
if !c.is_one() {
factors.push((self.constant(c), 1));
}
factors
}
fn factor(&self) -> Vec<(Self, usize)> {
if self.is_zero() {
return vec![];
}
let c = self.content();
let stripped = self.map_coeff(
|coeff| {
let coeff = self.ring.div(coeff, &c);
debug_assert!(coeff.is_integer());
coeff.numerator()
},
Z,
);
let mut factors: Vec<_> = stripped
.factor()
.into_iter()
.map(|(ff, p)| (ff.map_coeff(|coeff| coeff.into(), Q), p))
.collect();
if !c.is_one() {
factors.push((self.constant(c), 1));
}
factors
}
fn is_irreducible(&self) -> bool {
let c = self.content();
let stripped = self.map_coeff(
|coeff| {
let coeff = self.ring.div(coeff, &c);
debug_assert!(coeff.is_integer());
coeff.numerator()
},
Z,
);
stripped.is_irreducible()
}
}
impl<E: PositiveExponent> Factorize
for MultivariatePolynomial<AlgebraicExtension<RationalField>, E, LexOrder>
{
fn square_free_factorization(&self) -> Vec<(Self, usize)> {
if self.is_zero() {
return vec![];
}
let c = self.lcoeff();
let stripped = self.clone().make_monic();
let mut factors = vec![];
let fs = stripped.factor_separable();
for f in fs {
let mut nf = f.square_free_factorization_0_char();
factors.append(&mut nf);
}
if factors.is_empty() || !self.ring.is_one(&c) {
factors.insert(0, (self.constant(c), 1));
}
factors
}
/// Perform Trager's algorithm for factorization.
fn factor(&self) -> Vec<(Self, usize)> {
if self.is_zero() {
return vec![];
}
let sf = self.square_free_factorization();
let mut constant = self.ring.one();
let mut full_factors = vec![];
for (f, p) in &sf {
if f.is_constant() {
self.ring
.mul_assign(&mut constant, self.ring.pow(&f.get_constant(), *p as u64));
continue;
}
let (v, s, g, n) = f.norm_impl();
let mut factors = n.factor();
factors.retain(|(f, _)| !f.is_constant());
if factors.len() == 1 {
full_factors.push((f.clone(), *p));
continue;
}
let mut g_f = g.to_number_field(&self.ring);
for (f, b) in factors {
debug!("Rational factor {}", f);
let alpha_poly = g.variable(&self.get_vars_ref()[v]).unwrap()
+ g.variable(&self.ring.poly().variables[0]).unwrap()
* &g.constant((s as u64).into());
let f = f.to_number_field(&self.ring);
let gcd = f.gcd(&g_f);
g_f = g_f / &gcd;
let g = MultivariatePolynomial::from_number_field(&gcd)
.replace_with_poly(v, &alpha_poly)
.to_number_field(&self.ring);
let lc = g.lcoeff();
self.ring
.mul_assign(&mut constant, &self.ring.pow(&lc, (b * p) as u64));
full_factors.push((g.mul_coeff(self.ring.inv(&lc)), b * p));
}
}
if !self.ring.is_one(&constant) || full_factors.is_empty() {
full_factors.push((self.constant(constant), 1));
}
full_factors
}
fn is_irreducible(&self) -> bool {
// TODO: improve
self.factor().len() == 1
}
}
impl<
UField: FiniteFieldWorkspace,
F: GaloisField<Base = FiniteField<UField>> + PolynomialGCD<E>,
E: PositiveExponent,
> Factorize for MultivariatePolynomial<F, E, LexOrder>
where
FiniteField<UField>: Field + FiniteFieldCore<UField> + PolynomialGCD<u16>,
<FiniteField<UField> as Set>::Element: Copy,
AlgebraicExtension<<F as GaloisField>::Base>: PolynomialGCD<E>,
{
fn square_free_factorization(&self) -> Vec<(Self, usize)> {
let c = self.lcoeff();
let stripped = self.clone().make_monic();
let mut factors = vec![];
let fs = stripped.factor_separable();
for f in fs {
let mut nf = f.square_free_factorization_bernardin();
factors.append(&mut nf);
}
if factors.is_empty() || !self.ring.is_one(&c) {
factors.push((self.constant(c), 1))
}
factors
}
fn factor(&self) -> Vec<(Self, usize)> {
if self.is_zero() {
return vec![];
}
let sf = self.square_free_factorization();
let mut factors = vec![];
let mut degrees = vec![0; self.nvars()];
for (f, p) in sf {
debug!("SFF {} {}", f, p);
let mut var_count = 0;
for v in 0..self.nvars() {
degrees[v] = f.degree(v).to_u32() as usize;
if degrees[v] > 0 {
var_count += 1;
}
}
match var_count {
0 => {
factors.push((f, p));
}
1 => {
for (d2, f2) in f.distinct_degree_factorization() {
debug!("DDF {} {}", f2, d2);
for f3 in f2.equal_degree_factorization(d2) {
debug!("EDF {} {}", f3, p);
factors.push((f3, p));
}
}
}
2 => {
let mut order: Vec<_> = degrees
.iter()
.enumerate()
.filter(|(_, d)| **d > 0)
.collect();
order.sort_by_key(|o| Reverse(o.1));
let order: Vec<_> = order.into_iter().map(|(v, _)| v).collect();
factors.extend(
f.bivariate_factorization(order[0], order[1])
.into_iter()
.map(|ff| (ff, p)),
)
}
_ => {
// TODO: find better order
let mut order: Vec<_> = degrees
.iter()
.enumerate()
.filter(|(_, d)| **d > 0)
.collect();
order.sort_by_key(|o| Reverse(o.1));
let mut order: Vec<_> = order.into_iter().map(|(v, _)| v).collect();
factors.extend(
f.multivariate_factorization(&mut order, 10, None)
.into_iter()
.map(|ff| (ff, p)),
)
}
}
}
factors
}
fn is_irreducible(&self) -> bool {
let mut sf = self.square_free_factorization();
if sf.len() > 1 {
return false;
}
let (f, p) = sf.pop().unwrap();
let mut degrees = vec![0; self.nvars()];
debug!("SFF {} {}", f, p);
let mut var_count = 0;
for v in 0..self.nvars() {
degrees[v] = f.degree(v).to_u32() as usize;
if degrees[v] > 0 {
var_count += 1;
}
}
match var_count {
0 => true,
1 => {
let mut d = f.distinct_degree_factorization();
if d.len() > 1 {
return false;
}
let (d2, f2) = d.pop().unwrap();
f2.equal_degree_factorization(d2).len() == 1
}
2 => {
let mut order: Vec<_> = degrees
.iter()
.enumerate()
.filter(|(_, d)| **d > 0)
.collect();
order.sort_by_key(|o| Reverse(o.1));
let order: Vec<_> = order.into_iter().map(|(v, _)| v).collect();
f.bivariate_factorization(order[0], order[1]).len() == 1
}
_ => {
// TODO: find better order
let mut order: Vec<_> = degrees
.iter()
.enumerate()
.filter(|(_, d)| **d > 0)
.collect();
order.sort_by_key(|o| Reverse(o.1));
let mut order: Vec<_> = order.into_iter().map(|(v, _)| v).collect();
f.multivariate_factorization(&mut order, 10, None).len() == 1
}
}
}
}
impl<
UField: FiniteFieldWorkspace,
F: GaloisField<Base = FiniteField<UField>> + PolynomialGCD<E>,
E: PositiveExponent,
> MultivariatePolynomial<F, E, LexOrder>
where
FiniteField<UField>: Field + FiniteFieldCore<UField> + PolynomialGCD<u16>,
<FiniteField<UField> as Set>::Element: Copy,
AlgebraicExtension<<F as GaloisField>::Base>: PolynomialGCD<E>,
{
/// Bernardin's algorithm for square free factorization.
fn square_free_factorization_bernardin(&self) -> Vec<(Self, usize)> {
if self.is_constant() {
if self.is_one() {
return vec![];
} else {
return vec![(self.clone(), 1)];
}
}
let mut f = self.clone();
let mut h = HashMap::default();
let mut hr;
for var in 0..self.nvars() {
if f.degree(var) > E::zero() {
(f, hr) = f.square_free_factorization_ff_yun(var);
for (part, pow) in hr {
h.entry(pow)
.and_modify(|f| {
*f = &*f * ∂
})
.or_insert(part);
}
}
}
// take the pth root
// the coefficients remain unchanged, since x^1/p = x
// since the derivative in every var is 0, all powers are divisible by p
let p = self.ring.characteristic().to_u64().unwrap() as usize;
let mut b = f.clone();
for es in b.exponents_iter_mut() {
for e in es {
if e.is_zero() {
continue;
}
if p < u32::MAX as usize {
debug_assert_eq!(e.to_u32() as usize % p, 0);
*e = *e / E::from_u32(p as u32);
} else {
// at the moment exponents are limited to 32-bits
// so only the case where e = 0 is supported
assert!(*e == E::zero());
}
}
}
let mut factors = vec![];
let sub_factors = b.square_free_factorization_bernardin();
for (mut k, n) in sub_factors {
for (powh, hi) in &mut h {
if *powh < p {
let g = k.gcd(hi);
if !g.is_constant() {
k = k / &g;
*hi = &*hi / &g;
factors.push((g, n * p + *powh));
}
}
}
if !k.is_constant() {
factors.push((k, n * p));
}
}
for (powh, hi) in h {
if !hi.is_constant() {
factors.push((hi, powh));
}
}
factors
}
/// A modified version of Yun's square free factorization algorithm.
fn square_free_factorization_ff_yun(&self, var: usize) -> (Self, Vec<(Self, usize)>) {
let b = self.derivative(var);
let mut c = self.gcd(&b);
let mut w = self / &c;
let mut v = &b / &c;
let mut factors = vec![];
let mut i = 1;
while !w.is_constant() && i < self.ring.characteristic().to_u64().unwrap() as usize {
let z = v - w.derivative(var);
let g = w.gcd(&z);
w = w / &g;
v = z / &g;
c = c / &w;
if !g.is_one() {
factors.push((g, i));
}
i += 1
}
(c, factors)
}
/// Perform distinct degree factorization on a monic, univariate and square-free polynomial.
pub fn distinct_degree_factorization(&self) -> Vec<(usize, Self)> {
let Some(var) = self.last_exponents().iter().position(|x| *x > E::zero()) else {
return vec![(0, self.clone())]; // constant polynomial
};
let mut e = self.last_exponents().to_vec();
e[var] = E::one();
let x = self.monomial(self.ring.one(), e);
let mut factors = vec![];
let mut h = x.clone();
let mut f = self.clone();
let mut i: usize = 0;
while !f.is_one() {
i += 1;
h = h.exp_mod_univariate(self.ring.size().unwrap(), &mut f);
let mut g = f.gcd(&(&h - &x));
if !g.is_one() {
f = f.quot_rem_univariate(&mut g).0;
factors.push((i, g));
}
if f.last_exponents()[var] < E::from_u32(2 * (i as u32 + 1)) {
// f cannot be split more
if !f.is_constant() {
factors.push((f.last_exponents()[var].to_u32() as usize, f));
}
break;
}
}
factors
}
/// Perform Cantor-Zassenhaus's probabilistic algorithm for
/// finding irreducible factors of degree `d`.
pub fn equal_degree_factorization(&self, d: usize) -> Vec<Self> {
let s = self.clone().make_monic();
let Some(var) = self.last_exponents().iter().position(|x| *x > E::zero()) else {
if d == 1 {
return vec![s];
} else {
panic!("Degree mismatch for {self}: {d}");
}
};
let n = self.degree(var).to_u32() as usize;
if n == d {
return vec![s];
}
let mut rng = rng();
let mut random_poly = self.zero_with_capacity(d);
let mut exp = vec![E::zero(); self.nvars()];
let characteristic = self.ring.characteristic();
// compute inverse
let mut s_rev = s.clone();
s_rev.reverse();
let s_inv = s_rev.inverse_univariate(var, E::from_u32(n as u32 + 1));
let factor = loop {
// generate a random non-constant polynomial
random_poly.clear();
for i in 0..n {
let r = self
.ring
.sample(&mut rng, (0, characteristic.to_i64().unwrap_or(i64::MAX)));
if !self.ring.is_zero(&r) {
exp[var] = E::from_u32(i as u32);
random_poly.append_monomial(r, &exp);
}
}
if random_poly.degree(var) == E::zero() {
continue;
}
*random_poly.coefficients.last_mut().unwrap() = self.ring.one();
let g = random_poly.gcd(&s);
if !g.is_one() {
break g;
}
let b = if self.ring.characteristic() == 2 {
let max = self.ring.get_extension_degree() as usize * d;
let mut b = random_poly.clone();
let mut vcur = b.clone();
for _ in 1..max {
vcur = (&vcur * &vcur).rem(&s);
b = b + vcur.clone();
}
b
} else {
// TODO: use Frobenius map and modular composition to prevent computing large exponent poly^(p^d)
let p = self.ring.size().unwrap();
random_poly.exp_mod_univariate_fast(
var,
&(&p.pow(d as u64) - &1i64.into()) / &2i64.into(),
&s,
&s_inv,
) - self.one()
};
if b.is_constant() {
continue;
}
let g = b.gcd(&s);
if !g.is_one() && g != s {
break g;
}
};
let mut factors = factor.equal_degree_factorization(d);
factors.extend((self / &factor).equal_degree_factorization(d));
factors
}
/// Perform distinct and equal degree factorization on a square-free univariate polynomial.
fn factor_distinct_equal_degree(&self) -> Vec<Self> {
let mut factors = vec![];
for (d2, f2) in self.distinct_degree_factorization() {
debug!("DDF {} {}", f2, d2);
for f3 in f2.equal_degree_factorization(d2) {
debug!("EDF {}", f3);
factors.push(f3);
}
}
factors
}
/// Bernardin's algorithm based on
/// "A new bivariate Hensel lifting algorithm for n factors"
/// by Garrett Paluck. The formulation of the algorithm in other sources contain serious errors.
// TODO: merge with an almost similar method for the integer case. A modification that needs
// to be made here is to make the lcoeff_y=0 monic
fn bivariate_hensel_lift_bernardin(
&self,
interpolation_var: usize,
lcoeff: &Self,
univariate_factors: &[Self],
iterations: usize,
) -> Vec<Self> {
let y_poly = self.to_univariate_polynomial_list(interpolation_var);
// add the leading coefficient as a first factor
let mut factors = vec![lcoeff.replace(interpolation_var, &self.ring.zero())];
factors.extend_from_slice(univariate_factors);
// extract coefficients in y
let mut u: Vec<_> = factors
.iter()
.map(|f| {
let mut dense = vec![self.zero(); iterations + 1];
dense[0] = f.clone();
dense
})
.collect();
// update the first polynomial as it may contain y, since it's lcoeff
let y_lcoeff = lcoeff.to_univariate_polynomial_list(interpolation_var);
for (p, e) in y_lcoeff {
u[0][e.to_u32() as usize] = p;
}
let mut p = u.clone();
let mut cur_p = p[0][0].clone();
for x in &mut p.iter_mut().skip(1) {
cur_p = cur_p * &x[0];
x[0] = cur_p.clone();
}
let delta = Self::diophantine_univariate(&mut factors, &self.one());
for k in 1..iterations {
// extract the coefficient required to compute the error in y^k
// computed using a convolution
p[0][k] = u[0][k].clone();
for i in 1..factors.len() {
for j in 0..k {
p[i][k] = &p[i][k] + &(&p[i - 1][k - j] * &u[i][j]);
}
}
// find the kth power of y in f
// since we compute the error per power of y, we cannot stop on a 0 error
let e = if let Some((v, _)) = y_poly.iter().find(|e| e.1.to_u32() as usize == k) {
v - &p.last().unwrap()[k]
} else {
-p.last().unwrap()[k].clone()
};
if e.is_zero() {
continue;
}
for ((dp, f), d) in u.iter_mut().zip(factors.iter_mut()).zip(&delta) {
dp[k] = &dp[k] + &(d * &e).quot_rem_univariate(f).1;
}
// update the coefficients with the new y^k contributions
// note that the lcoeff[k] contribution is not new
let mut t = self.zero();
for i in 1..factors.len() {
t = &u[i][0] * &t + &u[i][k] * &p[i - 1][0];
p[i][k] = &p[i][k] + &t;
}
}
// convert dense polynomials to multivariate polynomials
u.into_iter()
.map(|ts| {
let mut new_poly = self.zero_with_capacity(ts.len());
for (i, mut f) in ts.into_iter().enumerate() {
for x in f.exponents_iter_mut() {
x[interpolation_var] = E::from_u32(i as u32);
}
new_poly = new_poly + f;
}
new_poly
})
.collect()
}
/// Compute the bivariate factorization of a square-free polynomial.
fn bivariate_factorization(&self, main_var: usize, interpolation_var: usize) -> Vec<Self> {
assert!(main_var != interpolation_var);
if self.bivariate_irreducibility_test() {
return vec![self.clone()];
}
// check for problems arising from canceling terms in the derivative
let der = self.derivative(main_var);
if der.is_zero() {
return self.bivariate_factorization(interpolation_var, main_var);
}
let g = self.gcd(&der);
if !g.is_constant() {
let mut factors = g.bivariate_factorization(main_var, interpolation_var);
factors.extend((self / &g).bivariate_factorization(main_var, interpolation_var));
return factors;
}
let mut sample_point = self.ring.zero();
let mut uni_f = self.replace(interpolation_var, &sample_point);
let mut i = 0;
let mut rng = rng();
loop {
if self.ring.size() == Some(i.into()) {
let field = self
.ring
.upgrade(self.ring.get_extension_degree().to_u64().unwrap() as usize + 1);
debug!(
"Upgrading to Galois field with exponent {}",
field.get_extension_degree()
);
let s_l = self.map_coeff(|c| self.ring.upgrade_element(c, &field), field.clone());
let facs = s_l.bivariate_factorization(main_var, interpolation_var);
return facs
.into_iter()
.map(|f| f.map_coeff(|c| self.ring.downgrade_element(c), self.ring.clone()))
.collect();
}
if self.degree(main_var) == uni_f.degree(main_var)
&& uni_f.gcd(&uni_f.derivative(main_var)).is_constant()
{
break;
}
// TODO: sample simple points first
sample_point = self.ring.sample(&mut rng, (0, i));
uni_f = self.replace(interpolation_var, &sample_point);
i += 1;
}
let mut d = self.degree(interpolation_var).to_u32();
let shifted_poly = if !self.ring.is_zero(&sample_point) {
self.shift_var_cached(interpolation_var, &sample_point)
} else {
self.clone()
};
let fs = uni_f.factor_distinct_equal_degree();
let mut lcoeff = shifted_poly.lcoeff_last_varorder(&[main_var, interpolation_var]);
let mut lc_d = lcoeff.degree(interpolation_var).to_u32();
let iter = (d + lc_d + 1) as usize;
let mut factors =
shifted_poly.bivariate_hensel_lift_bernardin(interpolation_var, &lcoeff, &fs, iter);
factors.swap_remove(0); // remove the lcoeff
let mut rec_factors = vec![];
// factor recombination
let mut s = 1;
let mut rest = shifted_poly;
'len: while 2 * s <= factors.len() {
let mut fs = CombinationIterator::new(factors.len(), s);
while let Some(cs) = fs.next() {
// TODO: multiply in the leading coefficient here,
// then we can skip the Pade approximation and reduce the
// number of iterations in the Hensel lifting to d + 1, like in the integer case?
let mut g = rest.constant(rest.lcoeff());
for (i, f) in factors.iter().enumerate() {
if cs.contains(&i) {
g = (&g * f).mod_var(interpolation_var, E::from_u32(iter as u32 + 1));
}
}
let y_polys: Vec<_> = g
.to_univariate_polynomial_list(main_var)
.into_iter()
.map(|(x, _)| x)
.collect();
let mut g_lcoeff = Self::lcoeff_reconstruct(&y_polys, d, lc_d);
g = (&g * &g_lcoeff)
.mod_var(interpolation_var, E::from_u32(d + 1))
.make_monic();
let (h, r) = rest.quot_rem(&g, true);
if r.is_zero() {
rec_factors.push(g);
for i in cs.iter().rev() {
factors.remove(*i);
}
rest = h;
lcoeff = lcoeff.quot_rem_univariate(&mut g_lcoeff).0;
lc_d = lcoeff.degree(interpolation_var).to_u32();
d = rest.degree(interpolation_var).to_u32();
continue 'len;
}
}
s += 1;
}
rec_factors.push(rest);
if !self.ring.is_zero(&sample_point) {
for x in &mut rec_factors {
// shift the polynomial to y - sample
*x = x.shift_var_cached(interpolation_var, &self.ring.neg(&sample_point));
}
}
rec_factors
}
/// Reconstruct the leading coefficient using a Pade approximation with numerator degree `deg_n` and
/// denominator degree `deg_d`. The resulting denominator should be a factor of the leading coefficient.
fn lcoeff_reconstruct(coeffs: &[Self], deg_n: u32, deg_d: u32) -> Self {
let mut lcoeff = coeffs[0].constant(coeffs[0].ring.one());
for x in coeffs {
let d = x.rational_approximant_univariate(deg_n, deg_d).unwrap().1;
if !d.is_one() {
let g = d.gcd(&lcoeff);
lcoeff = lcoeff * &(d / &g);
}
}
lcoeff
}
/// Sort the bivariate factors based on their univariate image so that they are
/// aligned between the different vars.
fn canonical_sort(
biv_polys: &[Self],
replace_var: usize,
sample_points: &[(usize, <F as Set>::Element)],
) -> Vec<(Self, <F as Set>::Element, Self)> {
let mut univariate_factors = biv_polys
.iter()
.map(|f| {
let mut u = f.clone();
for (v, p) in sample_points {
if *v == replace_var {
u = u.replace(*v, p);
}
}
(f.clone(), u.lcoeff(), u.make_monic())
})
.collect::<Vec<_>>();
univariate_factors.sort_by(|(_, _, a), (_, _, b)| {
a.exponents
.cmp(&b.exponents)
.then(a.coefficients.internal_cmp(&b.coefficients))
});
univariate_factors
}
/// Precompute the leading coefficients of the polynomial factors, using an
/// adapted version of Kaltofen's algorithm that has modifications of Martin Lee and Stanislav Poslavsky.
fn lcoeff_precomputation(
&self,
bivariate_factors: &[Self],
sample_points: &[(usize, <F as Set>::Element)],
order: &[usize],
) -> Result<(Vec<Self>, Vec<Self>), usize> {
let lcoeff = self.univariate_lcoeff(order[0]);
let sqf = lcoeff.square_free_factorization();
let mut lcoeff_square_free = self.one();
for (f, _) in &sqf {
lcoeff_square_free = lcoeff_square_free * f;
}
let sorted_main_factors = Self::canonical_sort(bivariate_factors, order[1], sample_points);
let mut true_lcoeffs: Vec<_> = sorted_main_factors
.iter()
.map(|(_, u, _)| self.constant(u.clone()))
.collect();
let main_bivariate_factors: Vec<_> =
sorted_main_factors.into_iter().map(|(f, _, _)| f).collect();
let mut lcoeff_left = lcoeff.clone();
for f in &true_lcoeffs {
lcoeff_left = lcoeff_left / f;
}
// TODO: smarter ordering
for (i, &var) in order[1..].iter().enumerate() {
if lcoeff_left.is_one() {
break;
}
if lcoeff_left.degree(var).is_zero() {
continue;
}
// only construct factors that depend on var
let c = lcoeff_square_free.univariate_content(var);
// make sure that the content removal does not change the unit
let mut c_eval = c.clone();
for (v, p) in sample_points {
c_eval = c_eval.replace(*v, p);
}
let lcoeff_square_free_pp = &lcoeff_square_free / &c * &c_eval;
debug!("Content-free lcsqf {}", lcoeff_square_free_pp);
// check if the evaluated leading coefficient remains square free
let mut poly_eval = lcoeff_square_free_pp.clone();
for (v, p) in sample_points {
if *v != var {
poly_eval = poly_eval.replace(*v, p);
}
}
let sqf = poly_eval.square_free_factorization();
if sqf.len() != 1 || sqf[0].1 != 1 {
debug!("Polynomial is not square free: {}", poly_eval);
return Err(main_bivariate_factors.len());
}
let bivariate_factors = if var == order[1] {
main_bivariate_factors.to_vec()
} else {
let mut poly_eval = self.clone();
for (v, p) in sample_points {
if *v != var {
poly_eval = poly_eval.replace(*v, p);
}
}
if poly_eval.degree(order[0]) != self.degree(order[0])
|| poly_eval.degree(var) != self.degree(var)
|| poly_eval.univariate_lcoeff(order[0]).degree(var) != lcoeff.degree(var)
{
debug!("Bad sample for reconstructing lcoeff: degrees do not match");
return Err(main_bivariate_factors.len());
}
let bivariate_factors: Vec<_> =
poly_eval.factor().into_iter().map(|(f, _)| f).collect();
if bivariate_factors.len() != main_bivariate_factors.len() {
return Err(bivariate_factors.len().min(main_bivariate_factors.len()));
}
Self::canonical_sort(&bivariate_factors, var, sample_points)
.into_iter()
.map(|(f, _, _)| f)
.collect()
};
let square_free_lc_biv_factors: Vec<_> = bivariate_factors
.iter()
.map(|f| f.univariate_lcoeff(order[0]).square_free_factorization())
.collect();
let basis = Self::gcd_free_basis(
square_free_lc_biv_factors
.iter()
.flatten()
.map(|x| x.0.clone())
.filter(|x| !x.is_constant())
.collect(),
);
if basis.is_empty() {
continue;
}
let lifted = if basis.len() == 1 {
vec![lcoeff_square_free_pp.clone()]
} else {
let mut new_order = order.to_vec();
new_order.swap(1, i + 1);
new_order.remove(0);
lcoeff_square_free_pp.multivariate_hensel_lift_with_auto_lcoeff_fixing(
&basis,
sample_points,
&new_order,
)
};
for (l, fac) in true_lcoeffs.iter_mut().zip(&square_free_lc_biv_factors) {
let mut contrib = self.one();
for (full, b) in lifted.iter().zip(&basis) {
// check if a GCD-free basis element is a factor of the leading coefficient of this bivariate factor
if let Some((_, m)) = fac.iter().find(|(f, _)| f == b || f.try_div(b).is_some())
{
for _ in 0..*m {
contrib = &contrib * full;
}
}
}
let g = contrib.gcd(l);
let mut new = contrib / &g;
// make sure the new part keeps the desired image coeff intact
let mut b_lc_eval = new.clone();
for (v, p) in sample_points {
b_lc_eval = b_lc_eval.replace(*v, p);
}
new = new / &b_lc_eval;
*l = &*l * &new;
lcoeff_left = &lcoeff_left / &new;
}
}
if !lcoeff_left.is_one() {
panic!(
"Could not reconstruct leading coefficient of {self}: order={order:?}, samples={sample_points:?} Rest = {lcoeff_left}"
);
}
Ok((main_bivariate_factors, true_lcoeffs))
}
fn multivariate_hensel_lift_with_auto_lcoeff_fixing(
&self,
factors: &[Self],
sample_points: &[(usize, <F as Set>::Element)],
order: &[usize],
) -> Vec<Self> {
let lcoeff = self.univariate_lcoeff(order[0]);
if lcoeff.is_constant() {
// the factors should be properly normalized
let (mut uni, delta) =
Self::univariate_diophantine_field(factors, order, sample_points);
return self
.multivariate_hensel_lifting(
factors,
&mut uni,
&delta,
sample_points,
None,
order,
1,
)
.unwrap();
}
// repeat the leading coefficient for every factor so that the leading coefficient is known
let padded_lcoeffs = vec![lcoeff.clone(); factors.len()];
let mut self_adjusted = self.clone();
for _ in 1..factors.len() {
self_adjusted = &self_adjusted * &lcoeff;
}
// set the proper lc
let mut lc_var_eval = lcoeff.clone();
for (v, p) in sample_points {
if *v != order[0] {
lc_var_eval = lc_var_eval.replace(*v, p);
}
}
let adjusted_factors: Vec<_> = factors
.iter()
.map(|f| f.clone().make_monic() * &lc_var_eval)
.collect();
let (mut uni, delta) =
Self::univariate_diophantine_field(&adjusted_factors, order, sample_points);
self_adjusted
.multivariate_hensel_lifting(
&adjusted_factors,
&mut uni,
&delta,
sample_points,
Some(&padded_lcoeffs),
order,
1,
)
.unwrap()
.into_iter()
.map(|f| {
let c = f.univariate_content(order[0]);
f / &c
})
.collect()
}
fn univariate_diophantine_field(
factors: &[Self],
order: &[usize],
sample_points: &[(usize, <F as Set>::Element)],
) -> (Vec<Self>, Vec<Self>) {
// produce univariate factors and univariate delta
let mut univariate_factors = factors.to_vec();
for f in &mut univariate_factors {
for (v, s) in sample_points {
if order[0] != *v {
*f = f.replace(*v, s);
}
}
}
let univariate_deltas = Self::diophantine_univariate(
&mut univariate_factors,
&factors[0].constant(factors[0].ring.one()),
);
(univariate_factors, univariate_deltas)
}
/// Perform multivariate factorization on a square-free polynomial.
fn multivariate_factorization(
&self,
order: &mut [usize],
mut coefficient_upper_bound: u64,
max_bivariate_factors: Option<usize>,
) -> Vec<Self> {
if let Some(m) = max_bivariate_factors
&& m == 1
{
return vec![self.clone()];
}
// check for problems arising from canceling terms in the derivative
let der = self.derivative(order[0]);
if der.is_zero() {
let mut new_order = order.to_vec();
let v = new_order.remove(0);
new_order.push(v);
return self.multivariate_factorization(
&mut new_order,
coefficient_upper_bound,
max_bivariate_factors,
);
}
let g = self.gcd(&der);
if !g.is_constant() {
let mut factors =
g.multivariate_factorization(order, coefficient_upper_bound, max_bivariate_factors);
factors.extend((self / &g).multivariate_factorization(
order,
coefficient_upper_bound,
max_bivariate_factors,
));
return factors;
}
// select a suitable evaluation point
let mut sample_points: Vec<_> = order[1..].iter().map(|i| (*i, self.ring.zero())).collect();
let mut uni_f;
let mut biv_f;
let mut rng = rng();
let degree = self.degree(order[0]);
let uni_lcoeff = self.univariate_lcoeff(order[0]);
let mut content_fail_count = 0;
let mut sample_fail = Integer::zero();
'new_sample: loop {
sample_fail += &1.into();
if &sample_fail * &2.into() > self.ring.size().unwrap() {
// the field is too small, upgrade
let field = self
.ring
.upgrade(self.ring.get_extension_degree().to_u64().unwrap() as usize + 1);
debug!(
"Upgrading to Galois field with exponent {}",
field.get_extension_degree()
);
let s_l = self.map_coeff(|c| self.ring.upgrade_element(c, &field), field.clone());
let facs = s_l.multivariate_factorization(
order,
coefficient_upper_bound,
max_bivariate_factors,
);
return facs
.into_iter()
.map(|f| f.map_coeff(|c| self.ring.downgrade_element(c), self.ring.clone()))
.collect();
}
for s in &mut sample_points {
s.1 = self
.ring
.nth(rng.random_range(0..=coefficient_upper_bound).into());
}
biv_f = self.clone();
for ((v, s), rem_var) in sample_points[1..].iter().zip(&order[1..]).rev() {
biv_f = biv_f.replace(*v, s);
if biv_f.degree(*rem_var) != self.degree(*rem_var) {
coefficient_upper_bound += 10;
continue 'new_sample;
}
}
// requirement for leading coefficient precomputation
if biv_f.univariate_lcoeff(order[0]).degree(order[1]) != uni_lcoeff.degree(order[1]) {
debug!(
"Degree of x{} in leading coefficient of bivariate image is wrong",
order[1]
);
coefficient_upper_bound += 10;
continue 'new_sample;
}
let biv_df = biv_f.derivative(order[0]);
uni_f = biv_f.replace(sample_points[0].0, &sample_points[0].1);
let uni_df = uni_f.derivative(order[0]);
if degree == biv_f.degree(order[0])
&& degree == uni_f.degree(order[0])
&& biv_f.gcd(&biv_df).is_constant()
&& uni_f.gcd(&uni_df).is_constant()
{
if !biv_f.univariate_content(order[0]).is_constant() {
content_fail_count += 1;
debug!("Univariate content is not constant");
if content_fail_count == 4 {
// it is likely that we will always find content for this variable ordering, so change the
// second variable
// TODO: is this guaranteed to work or should we also change the first variable?
let sec_var = order[1];
order.copy_within(2..order.len(), 1);
order[order.len() - 1] = sec_var;
for ((vs, _), v) in sample_points.iter_mut().zip(&order[1..]) {
*vs = *v;
}
debug!("Changed the second variable to {}", order[1]);
content_fail_count = 0;
}
} else {
break;
}
}
coefficient_upper_bound += 10;
}
for (v, s) in &sample_points {
debug!("Sample point {}={}", v, self.ring.printer(s));
}
let bivariate_factors = biv_f.bivariate_factorization(order[0], order[1]);
if bivariate_factors.len() == 1 {
// the polynomial is irreducible
return vec![self.clone()];
}
if let Some(max) = max_bivariate_factors
&& bivariate_factors.len() > max
{
return self.multivariate_factorization(
order,
coefficient_upper_bound,
max_bivariate_factors,
);
}
let (sorted_biv_factors, true_lcoeffs) =
match self.lcoeff_precomputation(&bivariate_factors, &sample_points, order) {
Ok((sorted_biv_factors, true_lcoeffs)) => (sorted_biv_factors, true_lcoeffs),
Err(max_biv) => {
// the leading coefficient computation failed because the bivaraite factorization was wrong
// try again with other sample points and a better bound
return self.multivariate_factorization(
order,
coefficient_upper_bound + 10,
Some(max_biv),
);
}
};
for (b, l) in sorted_biv_factors.iter().zip(&true_lcoeffs) {
debug!("Bivariate factor {} with true lcoeff {}", b, l);
}
let sorted_biv_factors = self.impose_true_lcoeffs_on_factors(
&sorted_biv_factors,
&true_lcoeffs,
&sample_points,
order,
2,
);
let (mut uni, delta) = MultivariatePolynomial::univariate_diophantine_field(
&sorted_biv_factors,
order,
&sample_points,
);
let factorization = match self.multivariate_hensel_lifting(
&sorted_biv_factors,
&mut uni,
&delta,
&sample_points,
Some(&true_lcoeffs),
order,
2,
) {
Some(factorization) => factorization,
None => {
return self.multivariate_factorization(
order,
coefficient_upper_bound + 10,
max_bivariate_factors,
);
}
};
// test the factorization
let mut test = self.one();
for f in &factorization {
debug!("Factor = {}", f);
test = &test * f;
}
if &test == self {
factorization
} else {
debug!(
"No immediate factorization of {} for sample points {:?}",
self, sample_points
);
// the bivariate factorization has too many factors, try again with other sample points
self.multivariate_factorization(
order,
coefficient_upper_bound + 10,
Some(max_bivariate_factors.unwrap_or(bivariate_factors.len()) - 1),
)
}
}
}
impl<F: Field, E: PositiveExponent> MultivariatePolynomial<F, E, LexOrder> {
fn multivariate_diophantine(
univariate_deltas: &[Self],
univariate_factors: &mut [Self],
prods: &[Self],
error: &Self,
order: &[usize],
sample_points: &[(usize, F::Element)],
degrees: &[usize],
mod_vars: &[MultivariatePolynomial<F, E, LexOrder>],
) -> Vec<Self> {
if order.len() == 1 {
return univariate_deltas
.iter()
.zip(univariate_factors)
.map(|(d, f)| (d * error).quot_rem_univariate(f).1)
.collect();
}
let last_var = *order.last().unwrap();
let shift = &sample_points.iter().find(|s| s.0 == last_var).unwrap().1;
let prods_mod = prods
.iter()
.map(|f| f.replace(last_var, shift))
.collect::<Vec<_>>();
let error_mod = error.replace(last_var, shift);
let previous_mod_vars = &mod_vars[..mod_vars.len().saturating_sub(1)];
debug!("dioph e[x{}^0] = {}", last_var, error_mod);
let mut deltas = Self::multivariate_diophantine(
univariate_deltas,
univariate_factors,
&prods_mod,
&error_mod,
&order[..order.len() - 1],
sample_points,
°rees[..order.len() - 1],
&mod_vars[..order.len() - 1],
);
for d in &mut deltas {
for m in mod_vars {
*d = d.quot_rem(m, false).1;
}
}
let mut exp = vec![E::zero(); error.nvars()];
exp[last_var] = E::one();
let var_pow = error
.monomial(error.ring.one(), exp)
.shift_var(last_var, &error.ring.neg(shift));
let mut cur_exponent;
let mut next_exponent = var_pow.clone();
for j in 1..=*degrees.last().unwrap() {
cur_exponent = next_exponent.clone();
next_exponent = &next_exponent * &var_pow;
let mut e = error.clone();
for (d, p) in deltas.iter().zip(prods) {
debug!("delta {} p {}", d, p);
e = &e - &(d * p);
for m in previous_mod_vars {
e = e.quot_rem(m, false).1;
}
// TODO: mod with (x-shift)^(j+1)?
// then we cannot break on 0 error
}
debug!("dioph e at x{}^{} = {}", last_var, j, e);
if e.is_zero() {
break;
}
// Take the jth Taylor coefficient in the current variable. Lower
// powers may survive as representatives in the ambient polynomial
// ring, but they do not contribute to this coefficient.
let shifted_e = e.shift_var(last_var, shift);
let mut e_mod = error.zero();
for (coefficient, degree) in shifted_e.to_univariate_polynomial_list(last_var) {
if degree.to_u32() as usize == j {
e_mod = e_mod + coefficient;
}
}
debug!("dioph e[x{}^{}] = {}", last_var, j, e_mod);
if e_mod.is_zero() {
continue;
}
let mut new_deltas = Self::multivariate_diophantine(
univariate_deltas,
univariate_factors,
&prods_mod,
&e_mod,
&order[..order.len() - 1],
sample_points,
°rees[..order.len() - 1],
&mod_vars[..order.len() - 1],
);
for (d, nd) in deltas.iter_mut().zip(&mut new_deltas) {
debug!("dioph d[x{}^{}] = {}", last_var, j, nd);
// multiply (y-s)^j and keep the accumulated correction reduced
let nd = &*nd * &cur_exponent;
*d = &*d + &nd;
for m in mod_vars {
*d = d.quot_rem(m, false).1;
}
}
}
deltas
}
fn sparse_multivariate_diophantine_from_skeleton(
factors: &[Self],
prods: &[Self],
error: &Self,
skeletons: &[Self],
order: &[usize],
) -> Option<Vec<Self>> {
if prods.len() != skeletons.len() {
return None;
}
if let Some(deltas) = Self::sparse_multivariate_diophantine_two_factor_by_sampling(
factors, prods, error, skeletons, order,
) {
return Some(deltas);
}
if factors.len() > 2 {
if let Some(deltas) = Self::sparse_multivariate_diophantine_by_sampling(
factors, prods, error, skeletons, order,
) {
return Some(deltas);
}
}
None
}
fn evaluate_using_exponents_univariate_grouped(
&self,
exp_evals: &[F::Element],
main_var: usize,
out: &mut Self,
) {
out.clear();
let degree = self.degree(main_var).to_u32() as usize;
let mut coefficients = vec![self.ring.zero(); degree + 1];
for (term, eval) in self.into_iter().zip(exp_evals) {
let degree = term.exponents[main_var].to_u32() as usize;
self.ring
.add_mul_assign(&mut coefficients[degree], term.coefficient, eval);
}
let mut exponent = vec![E::zero(); self.nvars()];
for (degree, coefficient) in coefficients.into_iter().enumerate() {
if !self.ring.is_zero(&coefficient) {
exponent[main_var] = E::from_u32(degree as u32);
out.coefficients.push(coefficient);
out.exponents.extend_from_slice(&exponent);
exponent[main_var] = E::zero();
}
}
}
fn sparse_multivariate_diophantine_two_factor_by_sampling(
factors: &[Self],
prods: &[Self],
error: &Self,
skeletons: &[Self],
order: &[usize],
) -> Option<Vec<Self>> {
if factors.len() != 2 || prods.len() != 2 || skeletons.len() != 2 || order.len() < 2 {
return None;
}
let characteristic = error.ring.characteristic();
if !characteristic.is_zero() && !characteristic.is_prime(0) {
return None;
}
let sparse_factor = skeletons
.iter()
.enumerate()
.filter(|(_, skeleton)| !skeleton.is_zero())
.min_by_key(|(_, skeleton)| skeleton.nterms())
.map(|(factor_index, _)| factor_index)?;
let dense_factor = 1 - sparse_factor;
let skeleton = &skeletons[sparse_factor];
let main_var = order[0];
let sample_vars = &order[1..];
let mut groups: Vec<(E, Vec<Vec<E>>, Vec<F::Element>)> = Vec::new();
for exponent in skeleton.exponents.chunks(skeleton.nvars()) {
let degree = exponent[main_var];
if let Some((_, exponents, _)) = groups.iter_mut().find(|(d, _, _)| *d == degree) {
exponents.push(exponent.to_vec());
} else {
groups.push((degree, vec![exponent.to_vec()], Vec::new()));
}
}
let samples_needed = groups
.iter()
.map(|(_, exponents, _)| exponents.len())
.max()
.unwrap_or(0);
if samples_needed == 0 || skeleton.nterms() > 512 {
return None;
}
let mut rng = rng();
'sample_base: for _ in 0..SPARSE_MDP_SAMPLE_BASE_ATTEMPTS {
let base_points = Self::sparse_interpolation_base_points(error, sample_vars, &mut rng);
for (_, exponents, sample_generators) in &mut groups {
sample_generators.clear();
let mut seen = HashSet::default();
for exponent in exponents {
let generator =
Self::evaluate_monomial_exponent(&error.ring, exponent, &base_points);
if !Self::sparse_interpolation_generator_is_usable(
error,
&generator,
sample_generators,
) || !seen.insert(generator.clone())
{
continue 'sample_base;
}
sample_generators.push(generator);
}
}
let mut rhs: Vec<Vec<F::Element>> = groups
.iter()
.map(|_| Vec::with_capacity(samples_needed))
.collect();
let mut cache = Self::sample_cache(error, factors, prods);
let error_base = error.evaluate_exponents(&base_points, &mut cache);
let factors_base = factors
.iter()
.map(|f| f.evaluate_exponents(&base_points, &mut cache))
.collect::<Vec<_>>();
let mut error_current = Cow::Borrowed(&error_base);
let mut factors_current = factors_base.iter().map(Cow::Borrowed).collect::<Vec<_>>();
let mut error_image =
error.zero_with_capacity(error.degree(main_var).to_u32() as usize + 1);
let mut factor_images = factors
.iter()
.map(|f| f.zero_with_capacity(f.degree(main_var).to_u32() as usize + 1))
.collect::<Vec<_>>();
for sample_index in 0..samples_needed {
if sample_index > 0 {
for (current, base) in error_current.to_mut().iter_mut().zip(&error_base) {
error.ring.mul_assign(current, base);
}
for (current, base) in factors_current.iter_mut().zip(&factors_base) {
for (c, b) in current.to_mut().iter_mut().zip(base) {
error.ring.mul_assign(c, b);
}
}
}
error.evaluate_using_exponents_univariate_grouped(
&error_current,
main_var,
&mut error_image,
);
for ((factor, current), image) in
factors.iter().zip(&factors_current).zip(&mut factor_images)
{
factor.evaluate_using_exponents_univariate_grouped(current, main_var, image);
}
let Some(deltas_image) =
Self::try_univariate_diophantine(&mut factor_images, &error_image)
else {
continue 'sample_base;
};
for (group_index, (degree, _, _)) in groups.iter().enumerate() {
rhs[group_index].push(Self::univariate_coefficient(
&deltas_image[sparse_factor],
main_var,
*degree,
));
}
}
let mut sparse_delta = skeleton.zero();
for ((_, exponents, sample_generators), rhs) in groups.iter().zip(&rhs) {
let coefficients = error.solve_shifted_transposed_vandermonde(
sample_generators,
&rhs[..exponents.len()],
);
for (coefficient, exponent) in coefficients.into_iter().zip(exponents) {
if !error.ring.is_zero(&coefficient) {
sparse_delta.append_monomial(coefficient, exponent);
}
}
}
let residual = error - &(&sparse_delta * &prods[sparse_factor]);
if residual
.ring
.try_inv(&prods[dense_factor].lcoeff())
.is_none()
{
continue;
}
let (dense_delta, rem) = residual.quot_rem(&prods[dense_factor], false);
if !rem.is_zero() {
continue;
}
let mut deltas = vec![error.zero(), error.zero()];
deltas[sparse_factor] = sparse_delta;
deltas[dense_factor] = dense_delta;
let mut check = error.zero();
for (delta, prod) in deltas.iter().zip(prods) {
check = check + delta * prod;
}
if &check == error {
return Some(deltas);
}
}
None
}
fn sparse_multivariate_diophantine_by_sampling(
factors: &[Self],
prods: &[Self],
error: &Self,
skeletons: &[Self],
order: &[usize],
) -> Option<Vec<Self>> {
if factors.len() != skeletons.len() || factors.len() != prods.len() || order.len() < 2 {
return None;
}
let main_var = order[0];
let sample_vars = &order[1..];
let mut groups = Vec::with_capacity(skeletons.len());
let mut samples_needed = 0usize;
let mut total_terms = 0usize;
for skeleton in skeletons {
let mut factor_groups: Vec<(E, Vec<Vec<E>>, Vec<F::Element>)> = Vec::new();
for exponent in skeleton.exponents.chunks(skeleton.nvars()) {
total_terms += 1;
let degree = exponent[main_var];
if let Some((_, exponents, _)) =
factor_groups.iter_mut().find(|(d, _, _)| *d == degree)
{
exponents.push(exponent.to_vec());
} else {
factor_groups.push((degree, vec![exponent.to_vec()], Vec::new()));
}
}
for (_, exponents, _) in &factor_groups {
samples_needed = samples_needed.max(exponents.len());
}
groups.push(factor_groups);
}
if samples_needed == 0 || total_terms > 512 {
return None;
}
let mut rng = rng();
'sample_base: for _ in 0..SPARSE_MDP_SAMPLE_BASE_ATTEMPTS {
let base_points = Self::sparse_interpolation_base_points(error, sample_vars, &mut rng);
for factor_groups in &mut groups {
for (_, exponents, sample_generators) in factor_groups {
sample_generators.clear();
let mut seen = HashSet::default();
for exponent in exponents {
let generator =
Self::evaluate_monomial_exponent(&error.ring, exponent, &base_points);
if !Self::sparse_interpolation_generator_is_usable(
error,
&generator,
sample_generators,
) || !seen.insert(generator.clone())
{
continue 'sample_base;
}
sample_generators.push(generator);
}
}
}
let mut rhs: Vec<Vec<Vec<F::Element>>> = groups
.iter()
.map(|factor_groups| {
factor_groups
.iter()
.map(|_| Vec::with_capacity(samples_needed))
.collect()
})
.collect();
let mut cache = Self::sample_cache(error, factors, prods);
let error_base = error.evaluate_exponents(&base_points, &mut cache);
let factors_base = factors
.iter()
.map(|f| f.evaluate_exponents(&base_points, &mut cache))
.collect::<Vec<_>>();
let mut error_current = Cow::Borrowed(&error_base);
let mut factors_current = factors_base.iter().map(Cow::Borrowed).collect::<Vec<_>>();
let mut error_image =
error.zero_with_capacity(error.degree(main_var).to_u32() as usize + 1);
let mut factor_images = factors
.iter()
.map(|f| f.zero_with_capacity(f.degree(main_var).to_u32() as usize + 1))
.collect::<Vec<_>>();
for sample_index in 0..samples_needed {
if sample_index > 0 {
for (current, base) in error_current.to_mut().iter_mut().zip(&error_base) {
error.ring.mul_assign(current, base);
}
for (current, base) in factors_current.iter_mut().zip(&factors_base) {
for (c, b) in current.to_mut().iter_mut().zip(base) {
error.ring.mul_assign(c, b);
}
}
}
error.evaluate_using_exponents_univariate_grouped(
&error_current,
main_var,
&mut error_image,
);
for ((factor, current), image) in
factors.iter().zip(&factors_current).zip(&mut factor_images)
{
factor.evaluate_using_exponents_univariate_grouped(current, main_var, image);
}
let Some(deltas_image) =
Self::try_univariate_diophantine(&mut factor_images, &error_image)
else {
continue 'sample_base;
};
for (factor_index, factor_groups) in groups.iter().enumerate() {
for (group_index, (degree, _, _)) in factor_groups.iter().enumerate() {
rhs[factor_index][group_index].push(Self::univariate_coefficient(
&deltas_image[factor_index],
main_var,
*degree,
));
}
}
}
let mut deltas: Vec<_> = skeletons.iter().map(|s| s.zero()).collect();
for (factor_index, factor_groups) in groups.iter().enumerate() {
for (group_index, (_, exponents, sample_generators)) in
factor_groups.iter().enumerate()
{
let coefficients = error.solve_shifted_transposed_vandermonde(
sample_generators,
&rhs[factor_index][group_index][..exponents.len()],
);
for (coefficient, exponent) in coefficients.into_iter().zip(exponents) {
if !error.ring.is_zero(&coefficient) {
deltas[factor_index].append_monomial(coefficient, exponent);
}
}
}
}
let mut check = error.zero();
for (delta, prod) in deltas.iter().zip(prods) {
check = check + delta * prod;
}
if &check == error {
return Some(deltas);
}
}
None
}
fn sparse_interpolation_base_points(
poly: &Self,
sample_vars: &[usize],
rng: &mut impl rand::RngCore,
) -> Vec<(usize, F::Element)> {
sample_vars
.iter()
.map(|v| {
let mut value = poly.ring.sample(rng, (1, MAX_RNG_PREFACTOR as i64));
let mut attempts = 0;
while poly.ring.is_zero(&value) && attempts < 8 {
value = poly.ring.sample(rng, (1, MAX_RNG_PREFACTOR as i64));
attempts += 1;
}
(*v, value)
})
.collect()
}
fn sparse_interpolation_generator_is_usable(
poly: &Self,
generator: &F::Element,
previous_generators: &[F::Element],
) -> bool {
if poly.ring.try_inv(generator).is_none() {
return false;
}
for prev in previous_generators {
let diff = poly.ring.sub(generator, prev);
if poly.ring.try_inv(&diff).is_none() {
return false;
}
}
true
}
fn evaluate_monomial_exponent(
ring: &F,
exponent: &[E],
sample_points: &[(usize, F::Element)],
) -> F::Element {
let mut value = ring.one();
for (var, sample) in sample_points {
let e = exponent[*var].to_u32();
if e > 0 {
ring.mul_assign(&mut value, &ring.pow(sample, e as u64));
}
}
value
}
fn sample_cache(f: &Self, factors: &[Self], prods: &[Self]) -> Vec<Vec<F::Element>> {
let mut degrees = (0..f.nvars())
.map(|var| f.degree(var).to_u32() as usize)
.collect::<Vec<_>>();
for p in factors.iter().chain(prods) {
for (var, degree) in degrees.iter_mut().enumerate() {
*degree = (*degree).max(p.degree(var).to_u32() as usize);
}
}
degrees
.into_iter()
.map(|degree| vec![f.ring.zero(); (degree + 1).min(POW_CACHE_SIZE)])
.collect()
}
fn univariate_coefficient(poly: &Self, var: usize, degree: E) -> F::Element {
let mut exponent = vec![E::zero(); poly.nvars()];
exponent[var] = degree;
poly.coefficient(&exponent)
.unwrap_or_else(|| poly.ring.zero())
}
fn try_univariate_diophantine(factors: &mut [Self], rhs: &Self) -> Option<Vec<Self>> {
if factors
.iter()
.any(|factor| factor.ring.try_inv(&factor.lcoeff()).is_none())
{
return None;
}
let mut cur = factors.last()?.clone();
let mut products = vec![cur.clone()];
for factor in factors[1..].iter().rev().skip(1) {
cur = cur * factor;
products.push(cur.clone());
}
products.reverse();
let mut deltas = Vec::with_capacity(factors.len());
let mut cur_s = rhs.clone();
for (factor, product) in factors.iter_mut().zip(&mut products) {
let (g, s, t) = factor.eea_univariate(product);
if !g.is_one() {
return None;
}
let new_s = (t * &cur_s).quot_rem_univariate(factor).1;
deltas.push(new_s);
cur_s = (s * &cur_s).quot_rem_univariate(product).1;
}
deltas.push(cur_s);
Some(deltas)
}
fn impose_true_lcoeffs_on_factors(
&self,
factors: &[Self],
true_lcoeffs: &[Self],
sample_points: &[(usize, F::Element)],
order: &[usize],
current_index: usize,
) -> Vec<Self> {
let mut factors_with_true_lcoeff = Vec::with_capacity(factors.len());
for (factor, true_lcoeff) in factors.iter().zip(true_lcoeffs) {
let mut lcoeff = true_lcoeff.clone();
for &var in &order[current_index + 1..] {
if let Some((_, sample)) = sample_points.iter().find(|(v, _)| *v == var) {
lcoeff = lcoeff.replace(var, sample);
}
}
let mut coefficients = factor.to_univariate_polynomial_list(order[0]);
coefficients.last_mut().unwrap().0 = lcoeff;
let mut fixed_factor = self.zero();
let mut exp = vec![E::zero(); self.nvars()];
for (coefficient, degree) in coefficients {
exp[order[0]] = degree;
fixed_factor = fixed_factor + coefficient.mul_exp(&exp);
}
factors_with_true_lcoeff.push(fixed_factor);
}
factors_with_true_lcoeff
}
fn multivariate_hensel_lifting(
&self,
factors: &[Self],
univariate_factors: &mut [Self],
univariate_deltas: &[Self],
sample_points: &[(usize, F::Element)],
true_lcoeffs: Option<&[Self]>,
order: &[usize],
start_index: usize,
) -> Option<Vec<Self>> {
debug!("Hensel lift {} with order {:?}", self, order);
let mut degrees: Vec<_> = order
.iter()
.map(|v| self.degree(*v).to_u32() as usize)
.collect();
let mut reconstructed_factors = factors.to_vec();
for v in start_index..order.len() {
// Replace the leading coefficient in x0 before this lift step.
let mut factors_with_true_lcoeff = if let Some(true_lcoeffs) = true_lcoeffs {
self.impose_true_lcoeffs_on_factors(
&reconstructed_factors,
true_lcoeffs,
sample_points,
order,
v,
)
} else {
reconstructed_factors
};
let mut f = self.clone();
for j in v + 1..order.len() {
for s in sample_points {
if s.0 == order[j] {
f = f.replace(s.0, &s.1);
}
}
}
// shift the polynomial such that the evaluation point is now at 0
// so that we can use a convolution for fast error computation
let shift = &sample_points.iter().find(|s| s.0 == order[v]).unwrap().1;
f = f.shift_var(order[v], shift);
for f in &mut factors_with_true_lcoeff {
*f = f.shift_var(order[v], shift);
}
let mut factor_products = vec![];
for i in 0..factors_with_true_lcoeff.len() {
// we cannot compute prod_i = prod_k(f_k)/f_i as the leading coefficient may
// not have an inverse in the almost-field p^k
let mut tot = self.one();
for (j, b) in factors_with_true_lcoeff.iter().enumerate() {
if i != j {
tot = tot * b;
}
}
factor_products.push(tot);
}
reconstructed_factors = f.multivariate_hensel_step(
univariate_deltas,
univariate_factors,
sample_points,
&mut factors_with_true_lcoeff,
&factor_products,
&order[..=v],
&mut degrees[..=v],
)?;
for f in &mut reconstructed_factors {
*f = f.shift_var(order[v], &self.ring.neg(shift));
}
for f in &reconstructed_factors {
debug!("Reconstructed factor {}", f);
}
}
Some(reconstructed_factors)
}
fn multivariate_hensel_step(
&self,
univariate_deltas: &[Self],
univariate_factors: &mut [Self],
sample_points: &[(usize, F::Element)],
factors: &mut [Self],
prods: &[Self],
order: &[usize],
degrees: &mut [usize],
) -> Option<Vec<Self>> {
let last_var = *order.last().unwrap();
let last_degree = *degrees.last().unwrap();
let y_poly = self.to_univariate_polynomial_list(last_var);
// extract coefficients in last_var
let mut u: Vec<_> = factors
.iter()
.map(|f| {
let mut dense = vec![self.zero(); last_degree + 1];
for (p, e) in f.to_univariate_polynomial_list(last_var) {
dense[e.to_u32() as usize] = p;
}
dense
})
.collect();
// TODO: do entire initialization here?
// the loop below cannot be cut short anyway, so it's not extra work to do it all here
let mut p = u.clone();
let mut cur_p = p[0][0].clone();
for x in &mut p.iter_mut().skip(1) {
for j in &mut *x {
*j = &*j * &cur_p;
}
cur_p = x[0].clone();
}
let prod_mod = prods
.iter()
.map(|f| f.replace(last_var, &self.ring.zero()))
.collect::<Vec<_>>();
let factors_mod = factors
.iter()
.map(|f| f.replace(last_var, &self.ring.zero()))
.collect::<Vec<_>>();
debug!("in shift {}", self);
debug!("deg {:?}", degrees);
// create the polynomials (x_i-shift_i)^deg used for modding during Hensel lifting
let mut mod_vars = Vec::with_capacity(order.len() - 2);
let mut exp = vec![E::zero(); self.nvars()];
for r in order[1..order.len() - 1]
.iter()
.zip(°rees[1..order.len() - 1])
{
let shift = &sample_points.iter().find(|s| s.0 == *r.0).unwrap().1;
exp[*r.0] = E::one();
let var_pow = self
.monomial(self.ring.one(), exp.clone())
.shift_var(*r.0, &self.ring.neg(shift))
.pow(r.1 + 1);
exp[*r.0] = E::zero();
mod_vars.push(var_pow);
}
for k in 1..=last_degree {
// extract the coefficient required to compute the error in y^k
// computed using a convolution
for i in 1..factors.len() {
for j in 0..k {
p[i][k] = &p[i][k] + &(&p[i - 1][k - j] * &u[i][j]);
}
}
// find the kth power of y in f
// since we compute the error per power of y, we cannot stop on a 0 error
let e = if let Some((v, _)) = y_poly.iter().find(|e| e.1.to_u32() as usize == k) {
v - &p.last().unwrap()[k]
} else {
-p.last().unwrap()[k].clone()
};
debug!("hensel e[x{}^{}] = {}", last_var, k, e);
if e.is_zero() {
continue;
}
let skeletons: Vec<_> = u.iter().map(|ui| ui[k - 1].clone()).collect();
let new_delta = Self::sparse_multivariate_diophantine_from_skeleton(
&factors_mod,
&prod_mod,
&e,
&skeletons,
&order[..order.len() - 1],
)
.unwrap_or_else(|| {
Self::multivariate_diophantine(
univariate_deltas,
univariate_factors,
&prod_mod,
&e,
&order[..order.len() - 1],
sample_points,
°rees[..order.len() - 1],
&mod_vars,
)
});
// update the coefficients with the new y^k contributions
let mut t = self.zero();
for (i, (du, d)) in u.iter_mut().zip(&new_delta).enumerate() {
debug!("hensel d[x{}^{}] = {}", last_var, k, d);
du[k] = &du[k] + d;
if i > 0 {
t = &du[0] * &t + d * &p[i - 1][0];
} else {
t = &t + d;
}
p[i][k] = &p[i][k] + &t;
}
}
// convert dense polynomials to multivariate polynomials
let lifted = u
.into_iter()
.map(|ts| {
let mut new_poly = self.zero_with_capacity(ts.len());
for (i, mut f) in ts.into_iter().enumerate() {
for x in f.exponents_iter_mut() {
debug_assert_eq!(x[last_var], E::zero());
x[last_var] = E::from_u32(i as u32);
}
new_poly = new_poly + f;
}
new_poly
})
.collect::<Vec<_>>();
Some(lifted)
}
}
impl<E: PositiveExponent> MultivariatePolynomial<IntegerRing, E, LexOrder> {
fn integer_factor_start_auto_decision(&self, order: &[usize]) -> (bool, f64) {
let bivariate_degree_space = order
.iter()
.take(2)
.map(|&var| self.degree(var).to_u32() as f64 + 1.0)
.product::<f64>();
let bivariate_box_density = if bivariate_degree_space == 0.0 {
0.0
} else {
self.nterms() as f64 / bivariate_degree_space
};
let use_bivariate =
bivariate_box_density <= INTEGER_FACTOR_BIVARIATE_SPARSE_BOX_DENSITY_THRESHOLD;
(use_bivariate, bivariate_box_density)
}
/// Hensel lift a solution of `self = u * w mod p` to `self = u * w mod max_p`
/// where `max_p` is a power of `p`.
///
/// If the lifting is successful, i.e. the error is 0 at some stage,
/// it will return `Ok((u,w))` where `u` and `w` are the true factors over
/// the integers. If a true factorization is not possible, it returns
/// `Err((u,w))` where `u` and `w` are monic.
pub fn hensel_lift<UField: FiniteFieldWorkspace>(
&self,
mut u: MultivariatePolynomial<FiniteField<UField>, E, LexOrder>,
mut w: MultivariatePolynomial<FiniteField<UField>, E, LexOrder>,
gamma: Option<Integer>,
max_p: &Integer,
) -> Result<(Self, Self), (Self, Self)>
where
FiniteField<UField>: Field + PolynomialGCD<E> + FiniteFieldCore<UField>,
Integer: ToFiniteField<UField>,
{
let lcoeff = self.lcoeff(); // lcoeff % p != 0
let mut gamma = gamma.unwrap_or(lcoeff.clone());
let lcoeff_p = lcoeff.to_finite_field(&u.ring);
let gamma_p = gamma.to_finite_field(&u.ring);
let field = u.ring.clone();
let p = field.get_prime().to_integer();
let a = self.clone().mul_coeff(gamma.clone());
u = u.make_monic().mul_coeff(gamma_p.clone());
w = w.make_monic().mul_coeff(lcoeff_p.clone());
let (_, s, t) = u.eea_univariate(&w);
debug_assert!((&s * &u + &t * &w).is_one());
let mut u_i = u.map_coeff(|c| field.to_symmetric_integer(c), Z);
let mut w_i = w.map_coeff(|c| field.to_symmetric_integer(c), Z);
// only replace the leading coefficient
*u_i.coefficients.last_mut().unwrap() = gamma.clone();
*w_i.coefficients.last_mut().unwrap() = lcoeff;
let mut e = &a - &(&u_i * &w_i);
let mut m = p.clone();
while !e.is_zero() && &m <= max_p {
let e_p = e.map_coeff(|c| (c / &m).to_finite_field(&field), field.clone());
let (q, r) = (&e_p * &s).quot_rem_univariate(&mut w);
let tau = &e_p * &t + q * &u;
u_i = u_i
+ tau
.map_coeff(|c| field.to_symmetric_integer(c), Z)
.mul_coeff(m.clone());
w_i = w_i
+ r.map_coeff(|c| field.to_symmetric_integer(c), Z)
.mul_coeff(m.clone());
e = &a - &(&u_i * &w_i);
m = &m * &p;
}
if e.is_zero() {
let content = u_i.content();
if !content.is_one() {
u_i = u_i.div_coeff(&content);
gamma = &gamma / &content;
}
if !gamma.is_one() {
w_i = w_i.div_coeff(&gamma); // true division is possible in this case
}
Ok((u_i, w_i))
} else {
if !u_i.lcoeff().is_one() {
let inv = u_i.lcoeff().mod_inverse(&m);
u_i = u_i.map_coeff(|c| (c * &inv).symmetric_mod(&m), Z);
}
if !w_i.lcoeff().is_one() {
let inv = w_i.lcoeff().mod_inverse(&m);
w_i = w_i.map_coeff(|c| (c * &inv).symmetric_mod(&m), Z);
}
Err((u_i, w_i))
}
}
/// Lift multiple factors by creating a binary tree and lifting each product.
fn multi_factor_hensel_lift(
&self,
hs: &[MultivariatePolynomial<Zp, E, LexOrder>],
max_p: &Integer,
) -> Vec<Self> {
if hs.len() == 1 {
if self.lcoeff().is_one() {
return vec![self.clone()];
} else {
let inv = self.lcoeff().mod_inverse(max_p);
let r = self.map_coeff(|c| (c * &inv).symmetric_mod(max_p), Z);
return vec![r];
}
}
let (gs, hs) = hs.split_at(hs.len() / 2);
let mut g = gs[0].one();
for x in gs {
g = g * x;
}
let mut h = hs[0].one();
for x in hs {
h = h * x;
}
let (g_i, h_i) = self.hensel_lift(g, h, None, max_p).unwrap_or_else(|e| e);
let mut factors = g_i.multi_factor_hensel_lift(gs, max_p);
factors.extend(h_i.multi_factor_hensel_lift(hs, max_p));
factors
}
/// Factor a square-free univariate polynomial over the integers by Hensel lifting factors computed over
/// a finite field image of the polynomial.
fn factor_reconstruct(&self) -> Vec<Self> {
let Some(var) = self.last_exponents().iter().position(|x| *x > E::zero()) else {
return vec![self.clone()]; // constant polynomial
};
let d = self.degree(var).to_u32();
if d == 1 {
return vec![self.clone()];
}
// select a suitable prime
// we try small primes first as the distinct and equal degree algorithms
// scale as log(p)
let mut field;
let mut f_p;
let mut pi = PrimeIteratorU64::new(101);
loop {
let p = pi.next().unwrap();
if p > u32::MAX as u64 {
panic!("Ran out of primes during factorization of {self}");
}
let p = p as u32;
if (&self.lcoeff() % &Integer::Single(p as i64)).is_zero() {
continue;
}
field = Zp::new(p);
f_p = self.map_coeff(|f| f.to_finite_field(&field), field.clone());
let df_p = f_p.derivative(var);
// check is f_p remains square-free
if f_p.gcd(&df_p).is_one() {
break;
}
}
let hs: Vec<_> = f_p.factor_distinct_equal_degree();
if hs.len() == 1 {
// the polynomial is irreducible
return vec![self.clone()];
}
let bound = self.coefficient_bound();
let p: Integer = (field.get_prime() as i64).into();
let mut max_p = p.clone();
while max_p < bound {
max_p = &max_p * &p;
}
let mut factors = self.multi_factor_hensel_lift(&hs, &max_p);
#[cfg(debug_assertions)]
for (h, h_p) in factors.iter().zip(&hs) {
let hh_p = h
.map_coeff(|c| c.to_finite_field(&field), field.clone())
.make_monic();
if &hh_p != h_p {
panic!("Mismatch of lifted factor: {hh_p} vs {h_p} in {self}");
}
}
let mut rec_factors = vec![];
// factor recombination
let mut s = 1;
let mut rest = self.clone();
'len: while 2 * s <= factors.len() {
let mut fs = CombinationIterator::new(factors.len(), s);
while let Some(cs) = fs.next() {
// check if the constant term matches
if rest.exponents[..rest.nvars()]
.iter()
.all(|e| *e == E::zero())
{
let mut g1 = rest.lcoeff();
let mut h1 = rest.lcoeff();
for (i, f) in factors.iter().enumerate() {
if f.exponents[..rest.nvars()].iter().all(|x| *x == E::zero()) {
if cs.contains(&i) {
g1 = (&g1 * &f.coefficients[0]).symmetric_mod(&max_p);
} else {
h1 = (&h1 * &f.coefficients[0]).symmetric_mod(&max_p);
}
}
}
// TODO: improve check
// for monic factors we can do &g1 * &h1 != &rest.lcoeff() * &rest.coefficients[0]
if (&g1 * &h1).abs() > bound {
continue;
}
}
let mut g = rest.constant(rest.lcoeff());
for (i, f) in factors.iter().enumerate() {
if cs.contains(&i) {
g = (&g * f).map_coeff(|i| i.clone().symmetric_mod(&max_p), Z);
}
}
let c = g.content();
g = g.div_coeff(&c);
let (h, r) = rest.quot_rem(&g, true);
if r.is_zero() {
// should always happen happen when |g1|_1 * |h1|_1 <= bound
rec_factors.push(g);
for i in cs.iter().rev() {
factors.remove(*i);
}
let c = h.content();
rest = h.div_coeff(&c);
continue 'len;
}
}
s += 1;
}
rec_factors.push(rest);
rec_factors
}
/// Lift a solution of `poly ≡ lcoeff * univariate_factors mod y mod p^k`
/// to `mod y^iterations mod p^k`.
///
/// Univariate factors must be monic and `lcoeff_y=0` should be as well.
fn bivariate_hensel_lift_bernardin(
poly: &MultivariatePolynomial<FiniteField<Integer>, E, LexOrder>,
interpolation_var: usize,
lcoeff: &MultivariatePolynomial<FiniteField<Integer>, E, LexOrder>,
univariate_factors: &[MultivariatePolynomial<FiniteField<Integer>, E, LexOrder>],
iterations: usize,
p: u32,
k: usize,
) -> Vec<MultivariatePolynomial<FiniteField<Integer>, E, LexOrder>> {
let finite_field = Zp::new(p);
// add the leading coefficient as a first factor
let mut factors = vec![lcoeff.replace(interpolation_var, &poly.ring.zero())];
for f in univariate_factors {
factors.push(f.clone());
}
let delta = Self::lift_diophantine_univariate(
&mut factors,
&poly.constant(poly.ring.one()),
finite_field.get_prime(),
k,
);
let y_poly = poly.to_univariate_polynomial_list(interpolation_var);
// extract coefficients in y
let mut u: Vec<_> = factors
.iter()
.map(|f| {
let mut dense = vec![poly.zero(); iterations + 1];
dense[0] = f.clone();
dense
})
.collect();
// update the first polynomial as it may contain y, since it's lcoeff
let y_lcoeff = lcoeff.to_univariate_polynomial_list(interpolation_var);
for (p, e) in y_lcoeff {
u[0][e.to_u32() as usize] = p;
}
let mut p = u.clone();
let mut cur_p = p[0][0].clone();
for x in &mut p.iter_mut().skip(1) {
cur_p = cur_p * &x[0];
x[0] = cur_p.clone();
}
for k in 1..iterations {
// extract the coefficient required to compute the error in y^k
// computed using a convolution
p[0][k] = u[0][k].clone();
for i in 1..factors.len() {
for j in 0..k {
p[i][k] = &p[i][k] + &(&p[i - 1][k - j] * &u[i][j]);
}
}
// find the kth power of y in f
// since we compute the error per power of y, we cannot stop on a 0 error
let e = if let Some((v, _)) = y_poly.iter().find(|e| e.1.to_u32() as usize == k) {
v - &p.last().unwrap()[k]
} else {
-p.last().unwrap()[k].clone()
};
if e.is_zero() {
continue;
}
for ((dp, f), d) in u.iter_mut().zip(factors.iter()).zip(&delta) {
dp[k] = &dp[k] + &(d * &e).quot_rem_univariate_monic(f).1;
}
// update the coefficients with the new y^k contributions
// note that the lcoeff[k] contribution is not new
let mut t = poly.zero();
for i in 1..factors.len() {
t = &u[i][0] * &t + &u[i][k] * &p[i - 1][0];
p[i][k] = &p[i][k] + &t;
}
}
// convert dense polynomials to multivariate polynomials
u.into_iter()
.map(|ts| {
let mut new_poly = poly.zero_with_capacity(ts.len());
for (i, mut f) in ts.into_iter().enumerate() {
for x in f.exponents_iter_mut() {
x[interpolation_var] = E::from_u32(i as u32);
}
new_poly = new_poly + f;
}
new_poly
})
.collect()
}
/// Factor a square-free bivariate polynomial over the integers.
fn bivariate_factor_reconstruct(&self, main_var: usize, interpolation_var: usize) -> Vec<Self> {
if self.bivariate_irreducibility_test() {
return vec![self.clone()];
}
let d2 = self.degree(interpolation_var).to_u32();
// select a suitable evaluation point, as small as possible as to not change the coefficient bound
let mut sample_point;
let mut uni_f;
let mut i = 0u64;
loop {
sample_point = i.into();
uni_f = self.replace(interpolation_var, &sample_point);
if self.degree(main_var) == uni_f.degree(main_var)
&& uni_f.gcd(&uni_f.derivative(main_var)).is_constant()
{
break;
}
i += 1;
}
// factor the univariate polynomial
let mut uni_fs: Vec<_> = uni_f
.factor()
.into_iter()
.map(|(f, p)| {
debug_assert_eq!(p, 1);
f
})
.collect();
// strip potential content
uni_fs.retain_mut(|f| !f.is_constant());
// select a suitable prime
// we try small primes first as the distinct and equal degree algorithms
// scale as log(p)
let mut pi = PrimeIteratorU64::new(101);
let mut field;
'new_prime: loop {
let p = pi.next().unwrap();
if p > u32::MAX as u64 {
panic!("Ran out of primes during factorization of {self}");
}
let p = p as u32;
if (&uni_f.lcoeff() % &Integer::Single(p as i64)).is_zero() {
continue;
}
field = Zp::new(p);
// make sure the factors stay coprime
let fs_p: Vec<_> = uni_fs
.iter()
.map(|f| f.map_coeff(|c| c.to_finite_field(&field), field.clone()))
.collect();
for (j, f) in fs_p.iter().enumerate() {
for g in &fs_p[j + 1..] {
if !f.gcd(g).is_one() {
continue 'new_prime;
}
}
}
break;
}
let shifted_poly = if !sample_point.is_zero() {
self.shift_var(interpolation_var, &sample_point)
} else {
self.clone()
};
// TODO: if bound is less than u64, we may also use Zp64 for the computation
let bound = shifted_poly.coefficient_bound();
let p = field.get_prime().to_integer();
let mut max_p = p.clone();
let mut k = 1;
while &max_p * 2 < bound {
max_p = &max_p * &p;
k += 1;
}
let mod_field = FiniteField::<Integer>::new(max_p.clone());
// make all factors monic, this is possible since the lcoeff is invertible mod p^k
let uni_fs_mod: Vec<_> = uni_fs
.iter()
.map(|f| {
let f1 = f.map_coeff(|c| mod_field.to_element(c.clone()), mod_field.clone());
f1.make_monic()
})
.collect();
let mut f_mod =
shifted_poly.map_coeff(|c| mod_field.to_element(c.clone()), mod_field.clone());
// make sure the lcoeff is monic at y=0
let inv_coeff = mod_field.inv(&mod_field.to_element(uni_f.lcoeff().clone()));
let f_mod_monic = f_mod.clone().mul_coeff(inv_coeff);
let lcoeff_monic = f_mod_monic.lcoeff_last_varorder(&[main_var, interpolation_var]);
let mut factors = Self::bivariate_hensel_lift_bernardin(
&f_mod_monic,
interpolation_var,
&lcoeff_monic,
&uni_fs_mod,
(d2 + 1) as usize,
field.get_prime(),
k,
);
factors.swap_remove(0); // remove the lcoeff
let mut rec_factors = vec![];
// factor recombination
let mut s = 1;
let mut lcoeff = f_mod.lcoeff_last_varorder(&[main_var, interpolation_var]);
let mut rest = shifted_poly;
'len: while 2 * s <= factors.len() {
let mut fs = CombinationIterator::new(factors.len(), s);
while let Some(cs) = fs.next() {
let mut g = lcoeff.clone();
for (i, f) in factors.iter().enumerate() {
if cs.contains(&i) {
g = (&g * f).mod_var(interpolation_var, E::from_u32(d2 + 1));
}
}
// convert to integer
let mut g_int = g.map_coeff(|c| mod_field.to_symmetric_integer(c), Z);
let content = g_int.univariate_content(main_var);
g_int = &g_int / &content;
let (h, r) = rest.quot_rem(&g_int, true);
if r.is_zero() {
rec_factors.push(g_int);
for i in cs.iter().rev() {
factors.remove(*i);
}
rest = h;
f_mod = rest.map_coeff(|c| mod_field.to_element(c.clone()), mod_field.clone());
lcoeff = f_mod.lcoeff_last_varorder(&[main_var, interpolation_var]);
continue 'len;
}
}
s += 1;
}
rec_factors.push(rest);
if !sample_point.is_zero() {
for x in &mut rec_factors {
// shift the polynomial to y - sample
*x = x.shift_var(interpolation_var, &self.ring.neg(&sample_point));
}
}
rec_factors
}
/// Solve a Diophantine equation over the ring `Z_p^k` using Newton iteration.
/// All factors must be monic.
fn lift_diophantine_univariate(
factors: &mut [MultivariatePolynomial<FiniteField<Integer>, E, LexOrder>],
rhs: &MultivariatePolynomial<FiniteField<Integer>, E, LexOrder>,
p: u32,
k: usize,
) -> Vec<MultivariatePolynomial<FiniteField<Integer>, E, LexOrder>> {
let field = Zp::new(p);
let prime: Integer = (p as u64).into();
let mut f_p: Vec<_> = factors
.iter()
.map(|f| {
f.map_coeff(
|c| rhs.ring.to_symmetric_integer(c).to_finite_field(&field),
field.clone(),
)
})
.collect();
let rhs_p = rhs.map_coeff(
|c| rhs.ring.to_symmetric_integer(c).to_finite_field(&field),
field.clone(),
);
// TODO: recycle from finite field computation that must have happened earlier
let mut delta =
MultivariatePolynomial::<Zp, E, LexOrder>::diophantine_univariate(&mut f_p, &rhs_p);
let mut deltas: Vec<_> = delta
.iter()
.map(|s| {
s.map_coeff(
|c| field.to_symmetric_integer(c).to_finite_field(&rhs.ring),
rhs.ring.clone(),
)
})
.collect();
if k == 1 {
return deltas;
}
let mut tot = rhs.constant(rhs.ring.one());
for f in factors.iter() {
tot = &tot * f;
}
let pi = factors
.iter_mut()
.map(|f| tot.quot_rem_univariate(f).0)
.collect::<Vec<_>>();
let mut m = prime.clone();
for _ in 1..k {
// TODO: is breaking on e=0 safe?
let mut e = rhs.clone();
for (dd, pp) in deltas.iter_mut().zip(&pi) {
e = &e - &(&*dd * pp);
}
let e_m = e.map_coeff(
|c| (&rhs.ring.to_symmetric_integer(c) / &m).to_finite_field(&field),
field.clone(),
);
for ((p, d_m), d) in f_p.iter_mut().zip(&mut delta).zip(deltas.iter_mut()) {
let new_delta = (&e_m * &*d_m).quot_rem_univariate(p).1;
*d = &*d
+ &new_delta.map_coeff(
|c| (&field.to_symmetric_integer(c) * &m).to_finite_field(&rhs.ring),
rhs.ring.clone(),
);
}
m = &m * ′
}
deltas
}
/// Compute the Gelfond bound for the coefficients magnitude of every factor of this polynomial
fn coefficient_bound(&self) -> Integer {
let max_norm = self.coefficients.iter().map(|x| x.abs()).max().unwrap();
let mut bound = Integer::one();
let mut total_degree = 0;
let mut non_zero_vars = 0;
for v in 0..self.nvars() {
let d = self.degree(v).to_u32() as u64;
if d > 0 {
non_zero_vars += 1;
total_degree += d;
bound *= &Integer::from(d + 1);
}
}
// move the 2^n into the sqrt to prevent precision loss when converting the sqrt
// to an integer
bound = &bound * &Integer::Single(2).pow((total_degree * 2).saturating_sub(non_zero_vars));
let root_bound = match &bound {
Integer::Single(b) => Integer::Single((*b as f64).sqrt() as i64),
Integer::Double(_) | Integer::Large(_) => bound.root(2),
};
bound = &root_bound + &1i64.into();
&bound * &(&max_norm * &self.lcoeff().abs())
}
/// Sort the bivariate factors based on their univariate image so that they are
/// aligned between the different vars.
// TODO: merge with the implementation for finite fields as the implementation
// is almost identical
fn canonical_sort(
biv_polys: &[Self],
replace_var: usize,
sample_points: &[(usize, Integer)],
) -> Vec<(Self, Integer, Self)> {
let mut univariate_factors = biv_polys
.iter()
.map(|f| {
let mut u = f.clone();
for (v, p) in sample_points {
if *v == replace_var {
u = u.replace(*v, p);
}
}
// make sure the representative is unique
let mut uni = u.clone().make_primitive();
if uni.lcoeff().is_negative() {
uni = -uni;
}
(f.clone(), u.lcoeff(), uni)
})
.collect::<Vec<_>>();
univariate_factors.sort_by(|(_, _, a), (_, _, b)| {
a.exponents
.cmp(&b.exponents)
.then(a.coefficients.partial_cmp(&b.coefficients).unwrap())
});
univariate_factors
}
fn reconstruct_lcoeffs_from_univariate_sample(
&self,
lcoeff: &Self,
lcoeff_factorization: &[(Self, usize)],
univariate_factors: &[Self],
sample_points: &[(usize, Integer)],
univariate_content: &Integer,
) -> Option<Vec<Self>> {
let mut lcoeff_factors = Vec::new();
for (f, pow) in lcoeff_factorization {
if f.is_constant() {
continue;
}
let mut eval = f.clone();
for (v, p) in sample_points {
eval = eval.replace(*v, p);
}
if !eval.is_constant() {
return None;
}
let eval = eval.get_constant().abs();
if eval <= 1 {
return None;
}
lcoeff_factors.push((f, pow, eval));
}
for (i, (_, _, a)) in lcoeff_factors.iter().enumerate() {
for (_, _, b) in &lcoeff_factors[i + 1..] {
if !a.gcd(b).is_one() {
return None;
}
}
}
let mut true_lcoeffs = vec![self.one(); univariate_factors.len()];
let mut used = vec![0usize; lcoeff_factors.len()];
let mut residual_lcs: Vec<_> = univariate_factors.iter().map(|f| f.lcoeff()).collect();
for (i, residual) in residual_lcs.iter_mut().enumerate() {
for (j, (factor, multiplicity, eval)) in lcoeff_factors.iter().enumerate() {
while used[j] < **multiplicity && (&residual.abs() % eval).is_zero() {
true_lcoeffs[i] = &true_lcoeffs[i] * factor;
*residual = &*residual / eval;
used[j] += 1;
}
}
}
if used
.iter()
.zip(&lcoeff_factors)
.any(|(used, (_, multiplicity, _))| *used != **multiplicity)
{
return None;
}
for (lcoeff_i, factor) in true_lcoeffs.iter_mut().zip(univariate_factors) {
let mut eval = lcoeff_i.clone();
for (v, p) in sample_points {
eval = eval.replace(*v, p);
}
if !eval.is_constant() {
return None;
}
let eval = eval.get_constant();
if eval.is_zero() {
return None;
}
let univariate_lcoeff = factor.lcoeff();
if (&univariate_lcoeff % &eval) != 0 {
return None;
}
*lcoeff_i = lcoeff_i.clone().mul_coeff(&univariate_lcoeff / &eval);
}
let mut product = self.one();
for l in &true_lcoeffs {
product = &product * l;
}
if product.clone().mul_coeff(univariate_content.clone()) == *lcoeff {
Some(true_lcoeffs)
} else {
None
}
}
#[allow(dead_code)]
fn wang_lcoeff_precomputation(
&self,
bivariate_factors: &[Self],
sample_points: &[(usize, Integer)],
order: &[usize],
) -> Result<(Integer, Vec<Self>, Vec<Self>), usize> {
let lcoeff = self.univariate_lcoeff(order[0]);
let lcoeff_factorization = lcoeff.factor();
let sorted_biv_factors = Self::canonical_sort(bivariate_factors, order[1], sample_points)
.into_iter()
.map(|(f, _, _)| f)
.collect::<Vec<_>>();
let univariate_factor_images = sorted_biv_factors
.iter()
.map(|factor| {
let mut image = factor.clone();
for (v, p) in sample_points {
image = image.replace(*v, p);
}
image
})
.collect::<Vec<_>>();
let Some(mut true_lcoeffs) = self.reconstruct_lcoeffs_from_univariate_sample(
&lcoeff,
&lcoeff_factorization,
&univariate_factor_images,
sample_points,
&Integer::one(),
) else {
return Err(sorted_biv_factors.len());
};
let mut sorted_biv_factors = sorted_biv_factors;
let mut lcoeff_left = self.one();
for (f, b) in true_lcoeffs.iter_mut().zip(&mut sorted_biv_factors) {
let mut b_eval = b.clone();
for (v, p) in sample_points {
b_eval = b_eval.replace(*v, p);
}
let b_lc = b_eval.lcoeff();
let mut f_eval = f.clone();
for (v, p) in sample_points {
f_eval = f_eval.replace(*v, p);
}
let f_lc = f_eval.lcoeff();
let lcm = b_lc.lcm(&f_lc);
let b_cor = &lcm / &b_lc;
let f_cor = lcm / &f_lc;
*b = b.clone().mul_coeff(b_cor);
lcoeff_left = lcoeff_left.div_coeff(&f_cor);
*f = f.clone().mul_coeff(f_cor);
}
Ok((lcoeff_left.get_constant(), sorted_biv_factors, true_lcoeffs))
}
fn wang_lcoeff_precomputation_univariate(
&self,
univariate_factors: &[Self],
sample_points: &[(usize, Integer)],
order: &[usize],
univariate_content: &Integer,
) -> Result<Vec<Self>, usize> {
let lcoeff = self.univariate_lcoeff(order[0]);
let lcoeff_factorization = lcoeff.factor();
self.reconstruct_lcoeffs_from_univariate_sample(
&lcoeff,
&lcoeff_factorization,
univariate_factors,
sample_points,
univariate_content,
)
.ok_or(univariate_factors.len())
}
/// Precompute the leading coefficients of the polynomial factors, using an
/// adapted version of Kaltofen's algorithm that has modifications of Martin Lee and Stanislav Poslavsky.
///
// TODO: merge with the implementation for finite fields as the implementation
// is almost identical
#[allow(dead_code)]
fn lcoeff_precomputation(
&self,
bivariate_factors: &[Self],
sample_points: &[(usize, Integer)],
order: &[usize],
bound: Integer,
p: u32,
k: usize,
) -> Result<(Integer, Vec<Self>, Vec<Self>), usize> {
let lcoeff = self.univariate_lcoeff(order[0]);
let sqf = lcoeff.square_free_factorization();
let mut lcoeff_square_free = self.one();
for (f, _) in &sqf {
lcoeff_square_free = &lcoeff_square_free * f;
}
let sorted_main_factors = Self::canonical_sort(bivariate_factors, order[1], sample_points);
let mut true_lcoeffs: Vec<_> = bivariate_factors.iter().map(|_| self.one()).collect();
let mut lcoeff_left = lcoeff.clone();
let mut main_bivariate_factors: Vec<_> =
sorted_main_factors.into_iter().map(|(f, _, _)| f).collect();
// TODO: smarter ordering
for (i, &var) in order[1..].iter().enumerate() {
if lcoeff_left.is_one() {
break;
}
if lcoeff_left.degree(var).is_zero() {
continue;
}
// only construct factors that depend on var and remove integer content and unit
let c = lcoeff_square_free.univariate_content(var);
let mut lcoeff_square_free_pp = &lcoeff_square_free / &c;
// check if the evaluated leading coefficient remains square free
let mut poly_eval = lcoeff_square_free_pp.clone();
for (v, p) in sample_points {
if *v != var {
poly_eval = poly_eval.replace(*v, p);
}
}
if poly_eval.lcoeff().is_negative() {
lcoeff_square_free_pp = -lcoeff_square_free_pp;
poly_eval = -poly_eval;
}
debug!("Content-free lcsqf {}", lcoeff_square_free_pp);
let sqf = poly_eval.square_free_factorization();
if sqf.len() != 1 || sqf[0].1 != 1 {
debug!("Polynomial is not square free: {}", poly_eval);
return Err(main_bivariate_factors.len());
}
let bivariate_factors = if var == order[1] {
main_bivariate_factors.to_vec()
} else {
let mut poly_eval = self.clone();
for (v, p) in sample_points {
if *v != var {
poly_eval = poly_eval.replace(*v, p);
}
}
if poly_eval.degree(order[0]) != self.degree(order[0])
|| poly_eval.degree(var) != self.degree(var)
|| poly_eval.univariate_lcoeff(order[0]).degree(var) != lcoeff.degree(var)
{
debug!("Bad sample for reconstructing lcoeff: degrees do not match");
return Err(main_bivariate_factors.len());
}
let bivariate_factors: Vec<_> = poly_eval
.factor()
.into_iter()
.map(|(f, _)| f)
// remove spurious content caused by particular evaluation point
.filter(|f| !f.is_constant())
.collect();
if bivariate_factors.len() != main_bivariate_factors.len() {
return Err(bivariate_factors.len().min(main_bivariate_factors.len()));
}
Self::canonical_sort(&bivariate_factors, var, sample_points)
.into_iter()
.map(|(f, _, _)| f)
.collect()
};
let square_free_lc_biv_factors: Vec<_> = bivariate_factors
.iter()
.map(|f| {
let mut sff = f.univariate_lcoeff(order[0]).square_free_factorization();
// make sure every bivariate factor has positive lcoeff such that the product
// of the basis elements equals the evaluated lcoeff_square_free_pp
for (b, _) in &mut sff {
if b.lcoeff().is_negative() {
*b = -b.clone();
}
}
sff
})
.collect();
let basis = Self::gcd_free_basis(
square_free_lc_biv_factors
.iter()
.flatten()
.map(|x| x.0.clone())
.filter(|x| !x.is_constant())
.collect(),
);
if basis.is_empty() {
continue;
}
let lifted = if basis.len() == 1 {
vec![lcoeff_square_free_pp.clone()]
} else {
let mut new_order = order.to_vec();
new_order.swap(1, i + 1);
new_order.remove(0);
lcoeff_square_free_pp.multivariate_hensel_lift_with_auto_lcoeff_fixing(
&basis,
sample_points,
&new_order,
bound.clone(),
p,
k,
)
};
for (l, fac) in true_lcoeffs.iter_mut().zip(&square_free_lc_biv_factors) {
let mut contrib = self.one();
for (full, b) in lifted.iter().zip(&basis) {
// check if a GCD-free basis element is a factor of the leading coefficient of this bivariate factor
if let Some((_, m)) = fac.iter().find(|(f, _)| f == b || f.try_div(b).is_some())
{
for _ in 0..*m {
contrib = &contrib * full;
}
}
}
let g = contrib.gcd(l);
let new = (contrib / &g).make_primitive();
*l = (&*l * &new).make_primitive();
let (q, r) = lcoeff_left.quot_rem(&new, true);
if !r.is_zero() {
panic!(
"Problem with bivariate factor scaling in factorization of {self}: order={order:?}, samples={sample_points:?}"
);
}
lcoeff_left = q;
}
}
if !lcoeff_left.is_constant() {
panic!(
"Could not reconstruct leading coefficient of {self}: order={order:?}, samples={sample_points:?} Rest = {lcoeff_left}"
);
}
// rescale the leading coefficient factors to recover the missing content and sign
for (f, b) in true_lcoeffs.iter_mut().zip(&mut main_bivariate_factors) {
let mut b_eval = b.clone();
for (v, p) in sample_points {
b_eval = b_eval.replace(*v, p);
}
let b_lc = b_eval.lcoeff();
let mut f_eval = f.clone();
for (v, p) in sample_points {
f_eval = f_eval.replace(*v, p);
}
let f_lc = f_eval.lcoeff();
let lcm = b_lc.lcm(&f_lc);
let b_cor = &lcm / &b_lc;
let f_cor = lcm / &f_lc;
*b = b.clone().mul_coeff(b_cor);
lcoeff_left = lcoeff_left.div_coeff(&f_cor);
*f = f.clone().mul_coeff(f_cor);
}
Ok((
lcoeff_left.get_constant(),
main_bivariate_factors,
true_lcoeffs,
))
}
#[allow(dead_code)]
fn multivariate_hensel_lift_with_auto_lcoeff_fixing(
&self,
factors: &[Self],
sample_points: &[(usize, Integer)],
order: &[usize],
bound: Integer,
p: u32,
k: usize,
) -> Vec<Self> {
let modulus = FiniteField::<Integer>::new(bound);
let ff = self.map_coeff(|c| modulus.to_element(c.clone()), modulus.clone());
let factors_ff: Vec<_> = factors
.iter()
.map(|f| f.map_coeff(|c| modulus.to_element(c.clone()), modulus.clone()))
.collect();
let sample_points_ff: Vec<_> = sample_points
.iter()
.map(|(v, p)| (*v, modulus.to_element(p.clone())))
.collect();
let lcoeff = ff.univariate_lcoeff(order[0]);
if lcoeff.is_constant() {
// the factors should be properly normalized
let (mut uni, delta) = MultivariatePolynomial::get_univariate_factors_and_deltas(
&factors_ff,
order,
sample_points,
p,
k,
);
let h = ff
.multivariate_hensel_lifting(
&factors_ff,
&mut uni,
&delta,
&sample_points_ff,
None,
order,
1,
)
.unwrap();
return h
.into_iter()
.map(|f| f.map_coeff(|c| modulus.to_symmetric_integer(c), Z))
.collect();
}
// repeat the leading coefficient for every factor so that the leading coefficient is known
let padded_lcoeffs = vec![lcoeff.clone(); factors.len()];
let mut self_adjusted = ff;
for _ in 1..factors_ff.len() {
self_adjusted = &self_adjusted * &lcoeff;
}
// set the proper lc
let mut lc_var_eval = lcoeff.clone();
for (v, p) in sample_points {
if *v != order[0] {
lc_var_eval = lc_var_eval.replace(*v, &lc_var_eval.ring.to_element(p.clone()));
}
}
let adjusted_factors: Vec<_> = factors_ff
.into_iter()
.map(|f| f.make_monic() * &lc_var_eval)
.collect();
let (mut uni, delta) = MultivariatePolynomial::get_univariate_factors_and_deltas(
&adjusted_factors,
order,
sample_points,
p,
k,
);
let h = self_adjusted
.multivariate_hensel_lifting(
&adjusted_factors,
&mut uni,
&delta,
&sample_points_ff,
Some(&padded_lcoeffs),
order,
1,
)
.unwrap();
h.into_iter()
.map(|f| {
let f_i = f.map_coeff(|c| modulus.to_symmetric_integer(c), Z);
let c = f_i.univariate_content(order[0]);
f_i / &c
})
.collect()
}
fn lcoeff_sample_preserves_square_free_images(
&self,
order: &[usize],
lcoeff: &Self,
lcoeff_square_free: &Self,
sample_points: &[(usize, Integer)],
) -> bool {
for &var in &order[1..] {
if lcoeff.degree(var).is_zero() {
continue;
}
let content = lcoeff_square_free.univariate_content(var);
let lcoeff_square_free_pp = lcoeff_square_free / &content;
let target_degree = lcoeff_square_free_pp.degree(var);
let mut poly_eval = lcoeff_square_free_pp;
for (v, p) in sample_points {
if *v != var {
poly_eval = poly_eval.replace(*v, p);
}
}
if poly_eval.lcoeff().is_negative() {
poly_eval = -poly_eval;
}
if poly_eval.degree(var) != target_degree {
return false;
}
let derivative = poly_eval.derivative(var);
if !poly_eval.gcd(&derivative).is_constant() {
return false; // not square-free
}
}
true
}
fn lcoeff_sample_supports_wang_reconstruction(
lcoeff_factorization: &[(Self, usize)],
sample_points: &[(usize, Integer)],
) -> bool {
let mut images = Vec::new();
for (factor, _) in lcoeff_factorization {
if factor.is_constant() {
continue;
}
let mut image = factor.clone();
for (v, p) in sample_points {
image = image.replace(*v, p);
}
if !image.is_constant() {
return false;
}
let image = image.get_constant().abs();
if image <= 1 {
return false;
}
if images
.iter()
.any(|prev: &Integer| !prev.gcd(&image).is_one())
{
return false;
}
images.push(image);
}
true
}
#[allow(dead_code)]
fn find_sample(
&self,
order: &mut [usize],
mut coefficient_upper_bound: i64,
mut max_factors_num: Option<usize>,
) -> (Vec<Self>, Vec<(usize, Integer)>, i64, Self) {
debug!("Find sample for {} with order {:?}", self, order);
// select a suitable evaluation point, as small as possible as to not change the coefficient bound
let mut cur_sample_points: Vec<_> =
order[1..].iter().map(|i| (*i, Integer::zero())).collect();
let mut cur_uni_f;
let mut cur_biv_f;
let mut rng = rng();
let degree = self.degree(order[0]);
let mut bivariate_factors: Vec<_>;
let mut best: Option<(Integer, Vec<Self>, Vec<(usize, Integer)>, i64, Self)> = None;
let uni_lcoeff = self.univariate_lcoeff(order[0]);
let mut lcoeff_square_free = self.one();
for (f, _) in uni_lcoeff.square_free_factorization() {
lcoeff_square_free = &lcoeff_square_free * &f;
}
let mut content_try_count = 0;
let mut lcoeff_try_count = 0;
'new_sample: loop {
for s in &mut cur_sample_points {
s.1 = Integer::Single(rng.random_range(0..=coefficient_upper_bound));
debug!("Sample x{} {}", s.0, s.1);
}
cur_biv_f = self.clone();
for ((v, s), rem_var) in cur_sample_points[1..].iter().zip(&order[1..]).rev() {
cur_biv_f = cur_biv_f.replace(*v, s);
if cur_biv_f.degree(*rem_var) != self.degree(*rem_var) {
coefficient_upper_bound += 10;
continue 'new_sample;
}
}
// requirement for leading coefficient precomputation
if cur_biv_f.univariate_lcoeff(order[0]).degree(order[1]) != uni_lcoeff.degree(order[1])
{
debug!(
"Degree of x{} in leading coefficient of bivariate image is wrong",
order[1]
);
coefficient_upper_bound += 10;
continue 'new_sample;
}
if !self.lcoeff_sample_preserves_square_free_images(
order,
&uni_lcoeff,
&lcoeff_square_free,
&cur_sample_points,
) {
debug!(
"Bad sample for reconstructing lcoeff: square-free lcoeff image is not square-free"
);
lcoeff_try_count += 1;
if lcoeff_try_count == 10 {
coefficient_upper_bound += 10;
lcoeff_try_count = 0;
}
continue 'new_sample;
}
lcoeff_try_count = 0;
let biv_df = cur_biv_f.derivative(order[0]);
cur_uni_f = cur_biv_f.replace(cur_sample_points[0].0, &cur_sample_points[0].1);
let uni_df = cur_uni_f.derivative(order[0]);
if degree == cur_biv_f.degree(order[0])
&& degree == cur_uni_f.degree(order[0])
&& cur_biv_f.gcd(&biv_df).is_constant()
&& cur_uni_f.gcd(&uni_df).is_constant()
{
let c = cur_biv_f.univariate_content(order[0]);
if !c.is_constant() {
content_try_count += 1;
coefficient_upper_bound += 10;
debug!("Univariate content is not constant");
if content_try_count == 10 {
// it is likely that we will always find content for this variable ordering, so change the
// second variable
// TODO: is this guaranteed to work or should we also change the first variable?
let sec_var = order[1];
order.copy_within(2..order.len(), 1);
order[order.len() - 1] = sec_var;
for ((vs, _), v) in cur_sample_points.iter_mut().zip(&order[1..]) {
*vs = *v;
}
debug!("Changed the second variable to {}", order[1]);
content_try_count = 0;
}
continue;
}
bivariate_factors = cur_biv_f.factor().into_iter().map(|f| f.0).collect();
bivariate_factors.retain(|f| !f.is_constant());
if max_factors_num.is_none() {
max_factors_num = Some(bivariate_factors.len());
}
// A valid specialization can split true multivariate factors, so the
// bivariate factor count is an upper bound on the final factor count.
// Keep the smallest admissible count; a larger count is only useful
// when the caller supplied a looser bound from an earlier failure.
if bivariate_factors.len() <= max_factors_num.unwrap() {
if bivariate_factors.len() < max_factors_num.unwrap() {
max_factors_num = Some(bivariate_factors.len());
content_try_count = 0;
best = None;
}
if best.is_none() || c.get_constant().abs() < best.as_ref().unwrap().0 {
best = Some((
c.get_constant().abs(),
bivariate_factors.clone(),
cur_sample_points.clone(),
coefficient_upper_bound,
cur_uni_f.clone(),
));
}
content_try_count += 1;
// try a few times to lower the chance of a costly Hensel lift
// with a wrong number of factors
if content_try_count > 2 {
break;
}
} else {
debug!(
"Number of factors is too large: {} vs {}",
bivariate_factors.len(),
max_factors_num.unwrap_or(bivariate_factors.len())
);
}
}
coefficient_upper_bound += 10;
debug!("Growing bound {}", coefficient_upper_bound);
}
let (_, bivariate_factors, cur_sample_points, coefficient_upper_bound, cur_uni_f) =
best.unwrap();
(
bivariate_factors,
cur_sample_points,
coefficient_upper_bound,
cur_uni_f,
)
}
fn find_univariate_sample(
&self,
order: &mut [usize],
mut coefficient_upper_bound: i64,
mut max_factors_num: Option<usize>,
) -> Option<(Vec<Self>, Vec<(usize, Integer)>, i64, Self, Integer)> {
debug!("Find univariate sample for {} with order {:?}", self, order);
let mut cur_sample_points: Vec<_> =
order[1..].iter().map(|i| (*i, Integer::zero())).collect();
let degree = self.degree(order[0]);
let mut seed = 0x9e37_79b9_7f4a_7c15u64;
for (i, &v) in order.iter().enumerate() {
seed ^= ((v as u64) + 0x517c_c1b7_2722_0a95u64).rotate_left(((i * 13) % 64) as u32);
}
seed ^= (self.nterms() as u64).rotate_left(17);
seed ^= (degree.to_u32() as u64).rotate_left(31);
let mut rng = StdRng::seed_from_u64(seed);
let mut best: Option<(usize, Vec<Self>, Vec<(usize, Integer)>, i64, Self, Integer)> = None;
let uni_lcoeff = self.univariate_lcoeff(order[0]);
let lcoeff_factorization = uni_lcoeff.factor();
let mut accepted_try_count = 0usize;
for _ in 0..512 {
for s in &mut cur_sample_points {
s.1 = Integer::Single(rng.random_range(0..=coefficient_upper_bound));
debug!("Sample x{} {}", s.0, s.1);
}
if !Self::lcoeff_sample_supports_wang_reconstruction(
&lcoeff_factorization,
&cur_sample_points,
) {
coefficient_upper_bound += 10;
continue;
}
let mut cur_uni_f = self.clone();
for (v, s) in cur_sample_points.iter().rev() {
cur_uni_f = cur_uni_f.replace(*v, s);
}
if degree != cur_uni_f.degree(order[0]) {
coefficient_upper_bound += 10;
continue;
}
let univariate_content = cur_uni_f.univariate_content(order[0]);
if !univariate_content.is_constant() {
coefficient_upper_bound += 10;
continue;
}
let univariate_content = univariate_content.get_constant();
let primitive_uni_f = cur_uni_f.clone().div_coeff(&univariate_content);
let uni_df = primitive_uni_f.derivative(order[0]);
if !primitive_uni_f.gcd(&uni_df).is_constant() {
coefficient_upper_bound += 10;
continue;
}
let mut univariate_factors: Vec<_> =
primitive_uni_f.factor().into_iter().map(|f| f.0).collect();
univariate_factors.retain(|f| !f.is_constant());
if univariate_factors.is_empty() {
coefficient_upper_bound += 10;
continue;
}
let mut product = self.one();
for f in &univariate_factors {
product = &product * f;
}
if product != primitive_uni_f {
if -product == primitive_uni_f {
univariate_factors[0] =
univariate_factors[0].clone().mul_coeff(Integer::from(-1));
} else {
coefficient_upper_bound += 10;
continue;
}
}
if self
.reconstruct_lcoeffs_from_univariate_sample(
&uni_lcoeff,
&lcoeff_factorization,
&univariate_factors,
&cur_sample_points,
&univariate_content,
)
.is_none()
{
coefficient_upper_bound += 10;
continue;
}
if max_factors_num.is_none() {
max_factors_num = Some(univariate_factors.len());
}
if univariate_factors.len() <= max_factors_num.unwrap() {
if univariate_factors.len() < max_factors_num.unwrap() {
max_factors_num = Some(univariate_factors.len());
accepted_try_count = 0;
best = None;
}
if best.as_ref().map_or(true, |(factor_count, _, _, _, _, _)| {
univariate_factors.len() < *factor_count
}) {
best = Some((
univariate_factors.len(),
univariate_factors.clone(),
cur_sample_points.clone(),
coefficient_upper_bound,
cur_uni_f.clone(),
univariate_content.clone(),
));
}
accepted_try_count += 1;
if accepted_try_count > 2 {
break;
}
}
coefficient_upper_bound += 10;
}
let Some((
_,
univariate_factors,
cur_sample_points,
coefficient_upper_bound,
cur_uni_f,
univariate_content,
)) = best
else {
return None;
};
Some((
univariate_factors,
cur_sample_points,
coefficient_upper_bound,
cur_uni_f,
univariate_content,
))
}
fn impose_true_lcoeffs_on_integer_factors(
&self,
factors: &[Self],
true_lcoeffs: &[Self],
order: &[usize],
) -> Vec<Self> {
let mut factors_with_true_lcoeff = Vec::with_capacity(factors.len());
for (factor, true_lcoeff) in factors.iter().zip(true_lcoeffs) {
let mut coefficients = factor.to_univariate_polynomial_list(order[0]);
coefficients.last_mut().unwrap().0 = true_lcoeff.clone();
let mut fixed_factor = self.zero();
let mut exp = vec![E::zero(); self.nvars()];
for (coefficient, degree) in coefficients {
exp[order[0]] = degree;
fixed_factor = fixed_factor + coefficient.mul_exp(&exp);
}
factors_with_true_lcoeff.push(fixed_factor);
}
factors_with_true_lcoeff
}
fn sparse_coefficient_hensel_lift_mod_prime(
&self,
mut factorization: Vec<Self>,
true_lcoeffs: &[Self],
p: u32,
max_p: &Integer,
order: &[usize],
) -> Option<Vec<Self>> {
let field = Zp::new(p);
let p_int: Integer = (p as u64).into();
factorization =
self.impose_true_lcoeffs_on_integer_factors(&factorization, true_lcoeffs, order);
let factors_mod_p: Vec<_> = factorization
.iter()
.map(|f| f.map_coeff(|c| c.to_finite_field(&field), field.clone()))
.collect();
let skeletons: Vec<_> = factors_mod_p
.iter()
.map(|f| f.mod_var(order[0], f.degree(order[0])))
.collect();
let prods_mod_p: Vec<_> = (0..factors_mod_p.len())
.map(|i| {
let mut prod = factors_mod_p[i].one();
for (j, f) in factors_mod_p.iter().enumerate() {
if i != j {
prod = prod * f;
}
}
prod
})
.collect();
let mut product = self.one();
for f in &factorization {
product = &product * f;
}
let mut error = self - &product;
let mut m = p_int.clone();
while !error.is_zero() && &m <= max_p {
let mut error_mod_p = factors_mod_p[0].zero();
for term in &error {
if !(term.coefficient % &m).is_zero() {
return None;
}
let q = term.coefficient / &m;
error_mod_p.append_monomial(q.to_finite_field(&field), term.exponents);
}
let deltas = MultivariatePolynomial::sparse_multivariate_diophantine_from_skeleton(
&factors_mod_p,
&prods_mod_p,
&error_mod_p,
&skeletons,
order,
)?;
for (factor, delta) in factorization.iter_mut().zip(deltas) {
*factor = &*factor
+ &delta
.map_coeff(|c| field.to_symmetric_integer(c), Z)
.mul_coeff(m.clone());
}
factorization =
self.impose_true_lcoeffs_on_integer_factors(&factorization, true_lcoeffs, order);
product = self.one();
for f in &factorization {
product = &product * f;
}
error = self - &product;
m = &m * &p_int;
}
Some(factorization)
}
fn multivariate_factorization_bivariate_start(
&self,
order: &mut [usize],
mut coefficient_upper_bound: i64,
mut max_bivariate_factors: Option<usize>,
) -> Vec<Self> {
if let Some(m) = max_bivariate_factors
&& m == 1
{
return vec![self.clone()];
}
let (bivariate_factors, sample_points, uni_f) = loop {
let (bivariate_factors, sample_points, coeff_b, uni_f) = self.find_sample(
order,
coefficient_upper_bound.max(10),
max_bivariate_factors,
);
coefficient_upper_bound = coeff_b;
if bivariate_factors.len() == 1 {
return vec![self.clone()];
}
if let Some(max) = max_bivariate_factors {
if bivariate_factors.len() < max {
debug!(
"Updating bivariate factor bound to {}",
bivariate_factors.len()
);
max_bivariate_factors = Some(bivariate_factors.len());
}
} else {
debug!(
"Updating bivariate factor bound to {}",
bivariate_factors.len()
);
max_bivariate_factors = Some(bivariate_factors.len());
}
break (bivariate_factors, sample_points, uni_f);
};
let mut prime_iter = PrimeIteratorU64::new(1 << 31);
let mut field;
let mut p;
'new_prime: loop {
p = prime_iter
.next()
.expect("Ran out of primes during factorization");
if p > u32::MAX as u64 {
panic!("Ran out of primes during factorization of {self}");
}
if (&uni_f.lcoeff() % &p.into()).is_zero() {
continue;
}
field = Zp64::new(p);
let fs_p: Vec<_> = bivariate_factors
.iter()
.map(|f| f.map_coeff(|c| c.to_finite_field(&field), field.clone()))
.collect();
for (f_p, f_z) in fs_p.iter().zip(&bivariate_factors) {
if f_p.degree(order[0]) != f_z.degree(order[0])
|| f_p.degree(order[1]) != f_z.degree(order[1])
|| f_p.ring.try_inv(&f_p.lcoeff()).is_none()
{
continue 'new_prime;
}
}
for (j, f) in fs_p.iter().enumerate() {
for g in &fs_p[j + 1..] {
if !f.gcd(g).is_one() {
continue 'new_prime;
}
}
}
break;
}
let p_int = field.get_prime().to_integer();
let mut lcoeff_max_p = p_int.clone();
let mut k = 1;
let lcoeff_bound = self.coefficient_bound();
while &lcoeff_max_p * 2 < lcoeff_bound {
lcoeff_max_p = &lcoeff_max_p * &p_int;
k += 1;
}
let (leftover_lc, mut sorted_biv_factors, mut true_lcoeffs) = match self
.lcoeff_precomputation(
&bivariate_factors,
&sample_points,
order,
lcoeff_max_p,
p as u32,
k,
) {
Ok((leftover_lc, sorted_biv_factors, true_lcoeffs)) => {
(leftover_lc, sorted_biv_factors, true_lcoeffs)
}
Err(max_biv) => {
return self.multivariate_factorization_bivariate_start(
order,
coefficient_upper_bound + 10,
Some(max_biv),
);
}
};
let rescaled = if leftover_lc != 1 {
for (b, l) in sorted_biv_factors.iter_mut().zip(&mut true_lcoeffs) {
*b = b.clone().mul_coeff(leftover_lc.clone());
*l = l.clone().mul_coeff(leftover_lc.clone());
}
Cow::Owned(
self.clone()
.mul_coeff(leftover_lc.pow(sorted_biv_factors.len() as u64 - 1)),
)
} else {
Cow::Borrowed(self)
};
let bound = rescaled.coefficient_bound();
let mut max_p = p_int.clone();
while &max_p * 2 < bound {
max_p = &max_p * &p_int;
}
for (b, l) in sorted_biv_factors.iter().zip(&true_lcoeffs) {
debug!("Bivariate factor {} with true lcoeff {}", b, l);
}
let poly_p = rescaled.map_coeff(|c| c.to_finite_field(&field), field.clone());
let sample_points_p: Vec<_> = sample_points
.iter()
.map(|(v, p)| (*v, p.to_finite_field(&field)))
.collect();
let true_lcoeffs_p: Vec<_> = true_lcoeffs
.iter()
.map(|f| f.map_coeff(|c| c.to_finite_field(&field), field.clone()))
.collect();
let sorted_biv_factors_p: Vec<_> = sorted_biv_factors
.iter()
.map(|f| f.map_coeff(|c| c.to_finite_field(&field), field.clone()))
.collect();
let sorted_biv_factors_p = poly_p.impose_true_lcoeffs_on_factors(
&sorted_biv_factors_p,
&true_lcoeffs_p,
&sample_points_p,
order,
2,
);
let (mut uni, delta) = MultivariatePolynomial::univariate_diophantine_field(
&sorted_biv_factors_p,
order,
&sample_points_p,
);
let Some(factorization_p) = poly_p.multivariate_hensel_lifting(
&sorted_biv_factors_p,
&mut uni,
&delta,
&sample_points_p,
Some(&true_lcoeffs_p),
order,
2,
) else {
return self.multivariate_factorization_bivariate_start(
order,
coefficient_upper_bound + 10,
max_bivariate_factors,
);
};
let factorization_z: Vec<_> = factorization_p
.into_iter()
.map(|f| f.map_coeff(|c| field.to_symmetric_integer(c), Z))
.collect();
let Some(factorization_z) = rescaled.sparse_coefficient_hensel_lift_mod_prime(
factorization_z,
&true_lcoeffs,
p as u32,
&max_p,
order,
) else {
return self.multivariate_factorization_bivariate_start(
order,
coefficient_upper_bound + 10,
max_bivariate_factors,
);
};
let mut factorization: Vec<MultivariatePolynomial<IntegerRing, E>> = factorization_z
.into_iter()
.map(|f| f.make_primitive())
.collect();
let mut test = self.one();
for f in &factorization {
debug!("Factor = {}", f);
test = &test * f;
}
if self.lcoeff().is_negative() != test.lcoeff().is_negative() {
test = -test;
if let Some(neg_coeff) = factorization.iter_mut().find(|f| f.lcoeff().is_negative()) {
*neg_coeff = -neg_coeff.clone();
} else {
factorization[0] = factorization[0].clone().mul_coeff((-1).into());
}
}
if &test == self {
let mut negated_factors = 0usize;
for f in &mut factorization {
if f.lcoeff().is_negative() {
*f = -f.clone();
negated_factors += 1;
}
}
if negated_factors % 2 == 1 {
factorization[0] = -factorization[0].clone();
}
factorization
} else {
self.multivariate_factorization_bivariate_start(
order,
coefficient_upper_bound + 10,
max_bivariate_factors,
)
}
}
/// Perform multivariate factorization on a square-free polynomial.
fn multivariate_factorization(
&self,
order: &mut [usize],
mut coefficient_upper_bound: i64,
mut max_univariate_factors: Option<usize>,
) -> Vec<Self> {
if let Some(m) = max_univariate_factors
&& m == 1
{
return vec![self.clone()];
}
let integer_start_mode = integer_factor_start_mode();
let use_bivariate_start = match integer_start_mode {
IntegerFactorStart::Auto => {
let (use_bivariate, _) = self.integer_factor_start_auto_decision(order);
use_bivariate
}
IntegerFactorStart::Univariate => false,
IntegerFactorStart::Bivariate => true,
IntegerFactorStart::Disabled => return vec![self.clone()],
};
if use_bivariate_start {
return self.multivariate_factorization_bivariate_start(
order,
coefficient_upper_bound,
max_univariate_factors,
);
}
let (univariate_factors, sample_points, uni_f, univariate_content) = loop {
let Some((univariate_factors, sample_points, coeff_b, uni_f, univariate_content)) =
self.find_univariate_sample(
order,
coefficient_upper_bound.max(10),
max_univariate_factors,
)
else {
if let Some(max) = max_univariate_factors {
panic!(
"failed to find a univariate factorization with at most {max} factors for order {order:?}; input={self}"
);
} else {
coefficient_upper_bound += 10;
}
continue;
};
coefficient_upper_bound = coeff_b;
if univariate_factors.len() == 1 {
// the polynomial is irreducible
return vec![self.clone()];
}
if let Some(max) = max_univariate_factors {
if univariate_factors.len() < max {
debug!(
"Updating univariate factor bound to {}",
univariate_factors.len()
);
max_univariate_factors = Some(univariate_factors.len());
}
} else {
debug!(
"Updating univariate factor bound to {}",
univariate_factors.len()
);
max_univariate_factors = Some(univariate_factors.len());
}
break (univariate_factors, sample_points, uni_f, univariate_content);
};
let univariate_factor_count = univariate_factors.len();
// select a suitable prime
let mut prime_iter = PrimeIteratorU64::new(1 << 31);
let mut field;
let mut p;
'new_prime: loop {
p = prime_iter
.next()
.expect("Ran out of primes during factorization");
if (&uni_f.lcoeff() % &p.into()).is_zero() {
continue;
}
field = Zp64::new(p);
// make sure the univariate factors stay square-free and coprime modulo p
let fs_p: Vec<_> = univariate_factors
.iter()
.map(|f| f.map_coeff(|c| c.to_finite_field(&field), field.clone()))
.collect();
for (f_p, f_z) in fs_p.iter().zip(&univariate_factors) {
if f_p.degree(order[0]) != f_z.degree(order[0])
|| f_p.ring.try_inv(&f_p.lcoeff()).is_none()
{
continue 'new_prime;
}
}
for (j, f) in fs_p.iter().enumerate() {
for g in &fs_p[j + 1..] {
if !f.gcd(g).is_one() {
continue 'new_prime;
}
}
}
break;
}
let lc_divs = match self.wang_lcoeff_precomputation_univariate(
&univariate_factors,
&sample_points,
order,
&univariate_content,
) {
Ok(true_lcoeffs) => true_lcoeffs,
Err(max_uni) => {
return self.multivariate_factorization(
order,
coefficient_upper_bound + 10,
Some(max_uni),
);
}
};
let scaled_true_lcoeffs: Vec<_> = lc_divs
.iter()
.map(|l| l.clone().mul_coeff(univariate_content.clone()))
.collect();
let scaled_univariate_factors: Vec<_> = univariate_factors
.iter()
.map(|f| f.clone().mul_coeff(univariate_content.clone()))
.collect();
let mut scale_pow = Integer::one();
for _ in 1..univariate_factor_count {
scale_pow *= &univariate_content;
}
let scaled_self = self.clone().mul_coeff(scale_pow);
let bound = scaled_self.coefficient_bound();
let p_int = field.get_prime().to_integer();
let mut max_p = p_int.clone();
while &max_p * 2 < bound {
max_p = &max_p * &p_int;
}
for (u, l) in univariate_factors.iter().zip(&lc_divs) {
debug!("Univariate factor {} with true lcoeff {}", u, l);
}
let poly_p = scaled_self.map_coeff(|c| c.to_finite_field(&field), field.clone());
let sample_points_p: Vec<_> = sample_points
.iter()
.map(|(v, p)| (*v, p.to_finite_field(&field)))
.collect();
let true_lcoeffs_p: Vec<_> = scaled_true_lcoeffs
.iter()
.map(|f| f.map_coeff(|c| c.to_finite_field(&field), field.clone()))
.collect();
let univariate_factors_p: Vec<_> = scaled_univariate_factors
.iter()
.map(|f| f.map_coeff(|c| c.to_finite_field(&field), field.clone()))
.collect();
let mut uni = univariate_factors_p.clone();
let Some(delta) = MultivariatePolynomial::try_univariate_diophantine(
&mut uni,
&univariate_factors_p[0].constant(field.one()),
) else {
return self.multivariate_factorization(
order,
coefficient_upper_bound + 10,
max_univariate_factors,
);
};
let Some(factorization_p) = poly_p.multivariate_hensel_lifting(
&univariate_factors_p,
&mut uni,
&delta,
&sample_points_p,
Some(&true_lcoeffs_p),
order,
1,
) else {
return self.multivariate_factorization(
order,
coefficient_upper_bound + 10,
max_univariate_factors,
);
};
let factorization_z: Vec<_> = factorization_p
.into_iter()
.map(|f| f.map_coeff(|c| field.to_symmetric_integer(c), Z))
.collect();
let Some(factorization_z) = scaled_self.sparse_coefficient_hensel_lift_mod_prime(
factorization_z,
&scaled_true_lcoeffs,
p as u32,
&max_p,
order,
) else {
return self.multivariate_factorization(
order,
coefficient_upper_bound + 10,
max_univariate_factors,
);
};
let mut factorization: Vec<MultivariatePolynomial<IntegerRing, E>> = factorization_z
.into_iter()
.map(|f| f.make_primitive())
.collect();
// test the factorization
let mut test = self.one();
for f in &factorization {
debug!("Factor = {}", f);
test = &test * f;
}
if self.lcoeff().is_negative() != test.lcoeff().is_negative() {
test = -test;
if let Some(neg_coeff) = factorization.iter_mut().find(|f| f.lcoeff().is_negative()) {
*neg_coeff = -neg_coeff.clone();
} else {
factorization[0] = factorization[0].clone().mul_coeff((-1).into());
}
}
if &test == self {
let mut negated_factors = 0usize;
for f in &mut factorization {
if f.lcoeff().is_negative() {
*f = -f.clone();
negated_factors += 1;
}
}
if negated_factors % 2 == 1 {
factorization[0] = -factorization[0].clone();
}
return factorization;
} else {
debug!(
"No immediate factorization of {} for sample points {:?}, retrying with factor bound {:?}",
self, sample_points, max_univariate_factors
);
return self.multivariate_factorization(
order,
coefficient_upper_bound + 10,
max_univariate_factors,
);
}
}
}
impl<E: PositiveExponent> MultivariatePolynomial<FiniteField<Integer>, E, LexOrder> {
/// Compute a univariate diophantine equation in `Z_p^k` by Newton iteration.
fn get_univariate_factors_and_deltas(
factors: &[Self],
order: &[usize],
sample_points: &[(usize, Integer)],
p: u32,
k: usize,
) -> (Vec<Self>, Vec<Self>) {
// produce univariate factors and univariate delta
let mut univariate_factors = factors.to_vec();
for f in &mut univariate_factors {
for (v, s) in sample_points {
if order[0] != *v {
*f = f.replace(*v, &f.ring.nth(s.clone()));
}
}
}
let univariate_deltas = MultivariatePolynomial::lift_diophantine_univariate(
&mut univariate_factors,
&factors[0].constant(factors[0].ring.one()),
p,
k,
);
(univariate_factors, univariate_deltas)
}
}
#[cfg(test)]
mod test {
use std::sync::{Arc, Mutex, atomic::Ordering};
use crate::{
GLOBAL_SETTINGS,
atom::AtomCore,
domains::{
InternalOrdering,
algebraic_number::AlgebraicExtension,
finite_field::{Z2, Zp},
integer::{Integer, Z},
rational::Q,
},
parse,
poly::factor::Factorize,
symbol,
};
static GLOBAL_FACTOR_SETTINGS_LOCK: Mutex<()> = Mutex::new(());
struct FactorSettingsGuard {
use_univariate: bool,
use_bivariate: bool,
}
impl FactorSettingsGuard {
fn new() -> Self {
Self {
use_univariate: GLOBAL_SETTINGS
.use_univariate_factorization
.load(Ordering::Relaxed),
use_bivariate: GLOBAL_SETTINGS
.use_bivariate_factorization
.load(Ordering::Relaxed),
}
}
}
impl Drop for FactorSettingsGuard {
fn drop(&mut self) {
GLOBAL_SETTINGS
.use_univariate_factorization
.store(self.use_univariate, Ordering::Relaxed);
GLOBAL_SETTINGS
.use_bivariate_factorization
.store(self.use_bivariate, Ordering::Relaxed);
}
}
#[test]
fn integer_factorization_respects_global_start_settings() {
let _lock = GLOBAL_FACTOR_SETTINGS_LOCK.lock().unwrap();
let _guard = FactorSettingsGuard::new();
let input = "(1+v1+v2+v3)*(2+3*v1+5*v2+7*v3)";
let poly = parse!(input).to_polynomial::<_, u8>(&Z, None);
GLOBAL_SETTINGS
.use_univariate_factorization
.store(true, Ordering::Relaxed);
GLOBAL_SETTINGS
.use_bivariate_factorization
.store(false, Ordering::Relaxed);
assert_eq!(poly.factor().len(), 2);
GLOBAL_SETTINGS
.use_univariate_factorization
.store(false, Ordering::Relaxed);
GLOBAL_SETTINGS
.use_bivariate_factorization
.store(true, Ordering::Relaxed);
assert_eq!(poly.factor().len(), 2);
GLOBAL_SETTINGS
.use_univariate_factorization
.store(false, Ordering::Relaxed);
GLOBAL_SETTINGS
.use_bivariate_factorization
.store(false, Ordering::Relaxed);
assert_eq!(poly.factor(), vec![(poly, 1)]);
}
#[test]
fn factor_ff_square_free() {
let field = Zp::new(3);
let poly = parse!("(1+v1)*(1+v1^2)^2*(v1^4+1)^3").to_polynomial::<_, u8>(&field, None);
let res = [("1+v1^4", 3), ("1+v1^2", 2), ("1+v1", 1)];
let mut res = res
.iter()
.map(|(f, p)| {
(
parse!(f)
.expand()
.to_polynomial(&field, poly.variables.clone()),
*p,
)
})
.collect::<Vec<_>>();
res.sort_by(|a, b| a.0.internal_cmp(&b.0).then(a.1.cmp(&b.1)));
let mut r = poly.square_free_factorization();
r.sort_by(|a, b| a.0.internal_cmp(&b.0).then(a.1.cmp(&b.1)));
assert_eq!(r, res);
}
#[test]
fn factor_ff_bivariate() {
let field = Zp::new(997);
let poly = parse!("((v2+1)*v1^2+v1*v2+1)*((v2^2+2)*v1^2+v2+1)")
.to_polynomial::<_, u8>(&field, None);
let res = [("1+2*v1^2+v2+v2^2*v1^2", 1), ("1+v1^2+v2*v1+v2*v1^2", 1)];
let mut res = res
.iter()
.map(|(f, p)| {
(
parse!(f)
.expand()
.to_polynomial(&field, poly.variables.clone()),
*p,
)
})
.collect::<Vec<_>>();
res.sort_by(|a, b| a.0.internal_cmp(&b.0).then(a.1.cmp(&b.1)));
let mut r = poly.factor();
r.sort_by(|a, b| a.0.internal_cmp(&b.0).then(a.1.cmp(&b.1)));
assert_eq!(r, res);
}
#[test]
fn factor_square_free() {
let poly =
parse!("3*(2*v1^2+v2)(v1^3+v2)^2(1+4*v2)^2(1+v1)").to_polynomial::<_, u8>(&Z, None);
let res = [
("3", 1),
("1+4*v2", 2),
("1+v1", 1),
("v2+2*v1^2", 1),
("v2+v1^3", 2),
];
let mut res = res
.iter()
.map(|(f, p)| {
(
parse!(f).expand().to_polynomial(&Z, poly.variables.clone()),
*p,
)
})
.collect::<Vec<_>>();
res.sort_by(|a, b| a.0.internal_cmp(&b.0).then(a.1.cmp(&b.1)));
let mut r = poly.square_free_factorization();
r.sort_by(|a, b| a.0.internal_cmp(&b.0).then(a.1.cmp(&b.1)));
assert_eq!(r, res);
}
#[test]
fn factor_univariate_1() {
let poly = parse!("2*(4 + 3*v1)*(3 + 2*v1 + 3*v1^2)*(3 + 8*v1^2)*(4 + v1 + v1^16)")
.to_polynomial::<_, u8>(&Z, None);
let res = [
("2", 1),
("4+3*v1", 1),
("3+2*v1+3*v1^2", 1),
("3+8*v1^2", 1),
("4+v1+v1^16", 1),
];
let mut res = res
.iter()
.map(|(f, p)| {
(
parse!(f).expand().to_polynomial(&Z, poly.variables.clone()),
*p,
)
})
.collect::<Vec<_>>();
res.sort_by(|a, b| a.0.internal_cmp(&b.0).then(a.1.cmp(&b.1)));
let mut r = poly.factor();
r.sort_by(|a, b| a.0.internal_cmp(&b.0).then(a.1.cmp(&b.1)));
assert_eq!(r, res);
}
#[test]
fn factor_univariate_2() {
let poly =
parse!("(v1+1)(v1+2)(v1+3)^3(v1+4)(v1+5)(v1^2+6)(v1^3+7)(v1+8)^2(v1^4+9)(v1^5+v1+10)")
.to_polynomial::<_, u8>(&Z, None);
let res = [
("5+v1", 1),
("1+v1", 1),
("4+v1", 1),
("2+v1", 1),
("7+v1^3", 1),
("10+v1+v1^5", 1),
("6+v1^2", 1),
("9+v1^4", 1),
("8+v1", 2),
("3+v1", 3),
];
let mut res = res
.iter()
.map(|(f, p)| {
(
parse!(f).expand().to_polynomial(&Z, poly.variables.clone()),
*p,
)
})
.collect::<Vec<_>>();
res.sort_by(|a, b| a.0.internal_cmp(&b.0).then(a.1.cmp(&b.1)));
let mut r = poly.factor();
r.sort_by(|a, b| a.0.internal_cmp(&b.0).then(a.1.cmp(&b.1)));
assert_eq!(r, res);
}
#[test]
fn factor_bivariate() {
let input = "(v1^2+v2+v1+1)(3*v1+v2^2+4)*(6*v1*(v2+1)+v2+5)*(7*v1*v2+4)";
let poly = parse!(input).to_polynomial::<_, u8>(&Z, None);
let res = [
("(1+v2+v1+v1^2)", 1),
("(5+v2+6*v1+6*v1*v2)", 1),
("(4+v2^2+3*v1)", 1),
("(4+7*v1*v2)", 1),
];
let mut res = res
.iter()
.map(|(f, p)| {
(
parse!(f).expand().to_polynomial(&Z, poly.variables.clone()),
*p,
)
})
.collect::<Vec<_>>();
res.sort_by(|a, b| a.0.internal_cmp(&b.0).then(a.1.cmp(&b.1)));
let mut r = poly.factor();
r.sort_by(|a, b| a.0.internal_cmp(&b.0).then(a.1.cmp(&b.1)));
assert_eq!(r, res);
}
#[test]
fn factor_multivariate() {
let input = "(v1*(2+2*v2+2*v3)+1)*(v1*(4+v3^2)+v2+3)*(v1*(v4+v4^2+4+v2)+v4+5)";
let poly = parse!(input).to_polynomial::<_, u8>(&Z, None);
let res = [
("5+v4+4*v1+v1*v4+v1*v4^2+v1*v2", 1),
("1+2*v1+2*v1*v3+2*v1*v2 ", 1),
("3+v2+4*v1+v1*v3^2", 1),
];
let mut res = res
.iter()
.map(|(f, p)| {
(
parse!(f).expand().to_polynomial(&Z, poly.variables.clone()),
*p,
)
})
.collect::<Vec<_>>();
res.sort_by(|a, b| a.0.internal_cmp(&b.0).then(a.1.cmp(&b.1)));
let mut r = poly.factor();
r.sort_by(|a, b| a.0.internal_cmp(&b.0).then(a.1.cmp(&b.1)));
assert_eq!(r, res);
}
#[test]
fn wang_lcoeff_reconstruction_accepts_unlucky_content_sample() {
let vars = Some(Arc::new(vec![
symbol!("m").into(),
symbol!("u").into(),
symbol!("v").into(),
]));
let poly = parse!("(3*m-u)*(5*m-v)")
.expand()
.to_polynomial::<_, u8>(&Z, vars);
let sample_points = vec![(1, Integer::from(3)), (2, Integer::from(10))];
let univariate_factors = ["m-1", "m-2"]
.iter()
.map(|f| parse!(f).to_polynomial(&Z, poly.variables.clone()))
.collect::<Vec<_>>();
let lcoeff = poly.univariate_lcoeff(0);
let lcoeff_factorization = lcoeff.factor();
assert!(
poly.reconstruct_lcoeffs_from_univariate_sample(
&lcoeff,
&lcoeff_factorization,
&univariate_factors,
&sample_points,
&Integer::from(1),
)
.is_none()
);
let reconstructed = poly
.reconstruct_lcoeffs_from_univariate_sample(
&lcoeff,
&lcoeff_factorization,
&univariate_factors,
&sample_points,
&Integer::from(15),
)
.unwrap();
assert_eq!(reconstructed, vec![poly.one(), poly.one()]);
}
const WANG_RESIDUAL_CONTENT_FACTORS: [&str; 4] = [
"3*m-u-w-v",
"128*m^6-128*m^5*u-352*m^5*w-112*m^5*v+40*m^4*u^2+232*m^4*u*w+100*m^4*u*v+272*m^4*w^2+432*m^4*w*v+32*m^4*v^2-4*m^3*u^3-30*m^3*u^2*w-22*m^3*u^2*v-134*m^3*u*w^2-276*m^3*u*w*v-25*m^3*u*v^2-80*m^3*w^3-266*m^3*w^2*v-210*m^3*w*v^2-3*m^3*v^3-2*m^2*u^3*w+m^2*u^3*v+12*m^2*u^2*w^2+33*m^2*u^2*w*v+3*m^2*u^2*v^2+22*m^2*u*w^3+126*m^2*u*w^2*v+123*m^2*u*w*v^2+2*m^2*u*v^3+8*m^2*w^4+50*m^2*w^3*v+89*m^2*w^2*v^2+47*m^2*w*v^3+m*u^3*w*v-8*m*u^2*w^2*v-11*m*u^2*w*v^2-13*m*u*w^3*v-38*m*u*w^2*v^2-25*m*u*w*v^3-2*m*w^4*v-8*m*w^3*v^2-10*m*w^2*v^3-4*m*w*v^4+u^2*w^2*v^2+u^2*w*v^3+2*u*w^3*v^2+4*u*w^2*v^3+2*u*w*v^4",
"128*m^4-80*m^3*u-80*m^3*w-96*m^3*v+12*m^2*u^2+52*m^2*u*v+12*m^2*w^2+52*m^2*w*v+24*m^2*v^2+9*m*u^2*w-10*m*u^2*v+9*m*u*w^2+4*m*u*w*v-12*m*u*v^2-10*m*w^2*v-12*m*w*v^2-2*m*v^3+u^3*v-3*u^2*w*v+2*u^2*v^2-3*u*w^2*v-2*u*w*v^2+u*v^3+w^3*v+2*w^2*v^2+w*v^3",
"16*m^3-8*m^2*u-8*m^2*w-20*m^2*v+m*u^2+2*m*u*w+6*m*u*v+m*w^2+6*m*w*v+8*m*v^2-u*w*v-u*v^2-w*v^2-v^3",
];
fn wang_residual_content_input() -> String {
WANG_RESIDUAL_CONTENT_FACTORS
.iter()
.map(|f| format!("({f})"))
.collect::<Vec<_>>()
.join("*")
}
#[test]
fn wang_univariate_lcoeff_keeps_residual_content_separate() {
let vars = Some(Arc::new(vec![
symbol!("m").into(),
symbol!("u").into(),
symbol!("w").into(),
symbol!("v").into(),
]));
let poly = parse!(wang_residual_content_input())
.expand()
.to_polynomial::<_, u16>(&Z, vars);
let mut order = vec![0, 3, 2, 1];
let (univariate_factors, sample_points, _, _, univariate_content) =
poly.find_univariate_sample(&mut order, 10, None).unwrap();
let lcoeff_divs = poly
.wang_lcoeff_precomputation_univariate(
&univariate_factors,
&sample_points,
&order,
&univariate_content,
)
.unwrap();
let lcoeff_div_product = lcoeff_divs
.iter()
.fold(poly.one(), |product, lcoeff| &product * lcoeff);
let actual_lcoeff = poly.univariate_lcoeff(order[0]);
assert!(!univariate_content.is_one());
assert_ne!(lcoeff_div_product, actual_lcoeff);
assert_eq!(
lcoeff_div_product.mul_coeff(univariate_content),
actual_lcoeff
);
}
#[test]
fn factor_multivariate_wang_residual_content() {
let vars = Some(Arc::new(vec![
symbol!("m").into(),
symbol!("u").into(),
symbol!("w").into(),
symbol!("v").into(),
]));
let poly = parse!(wang_residual_content_input())
.expand()
.to_polynomial::<_, u16>(&Z, vars);
let mut res = WANG_RESIDUAL_CONTENT_FACTORS
.iter()
.map(|f| {
(
parse!(f).expand().to_polynomial(&Z, poly.variables.clone()),
1,
)
})
.collect::<Vec<_>>();
res.sort_by(|a, b| a.0.internal_cmp(&b.0).then(a.1.cmp(&b.1)));
let mut r = poly.factor();
r.sort_by(|a, b| a.0.internal_cmp(&b.0).then(a.1.cmp(&b.1)));
assert_eq!(r, res);
}
#[test]
fn factor_overall_minus() {
let poly = parse!("-v1*v3^2-v1*v2*v3^2").to_polynomial::<_, u8>(
&Z,
Some(Arc::new(vec![
symbol!("v1").into(),
symbol!("v2").into(),
symbol!("v3").into(),
])),
);
let res = [("-1", 1), ("v3", 2), ("1+v2", 1), ("v1", 1)];
let mut res = res
.iter()
.map(|(f, p)| {
(
parse!(f).expand().to_polynomial(&Z, poly.variables.clone()),
*p,
)
})
.collect::<Vec<_>>();
res.sort_by(|a, b| a.0.internal_cmp(&b.0).then(a.1.cmp(&b.1)));
let mut r = poly.factor();
r.sort_by(|a, b| a.0.internal_cmp(&b.0).then(a.1.cmp(&b.1)));
assert_eq!(r, res);
}
#[test]
fn factor_multivariate_2() {
let poly = parse!("v2^2*v3-v1*v2*v3+v1*v2*v3^2+v1*v2^2-v1^2*v3^2+v1^2*v2*v3")
.to_polynomial::<_, u8>(
&Z,
Some(Arc::new(vec![
symbol!("v1").into(),
symbol!("v2").into(),
symbol!("v3").into(),
])),
);
let res = [("v2+v1*v3", 1), ("v2*v3-v1*v3+v1*v2", 1)];
let mut res = res
.iter()
.map(|(f, p)| {
(
parse!(f).expand().to_polynomial(&Z, poly.variables.clone()),
*p,
)
})
.collect::<Vec<_>>();
res.sort_by(|a, b| a.0.internal_cmp(&b.0).then(a.1.cmp(&b.1)));
let mut r = poly.factor();
r.sort_by(|a, b| a.0.internal_cmp(&b.0).then(a.1.cmp(&b.1)));
assert_eq!(r, res);
}
#[test]
fn galois_upgrade() {
let a =
parse!("x^7(y^5+y^4+y^3+y^2)+x^5(y^3+y)+x^4(y^4+y)+x^3(y^2+y)+x^2y+x*y^2+x*y+x+y+1")
.to_polynomial::<_, u8>(&Z2, None);
assert_eq!(a.factor().len(), 2)
}
#[test]
fn algebraic_extension() {
let a = parse!("z^4+z^3+(2+a-a^2)z^2+(1+a^2-2a^3)z-2").to_polynomial::<_, u8>(&Q, None);
let f = parse!("a^4-3").to_polynomial::<_, u16>(&Q, None);
let f = AlgebraicExtension::new(f);
let mut factors = a.to_number_field(&f).factor();
let f1 = parse!("(1-a^2)+(1-a)*z+z^2")
.to_polynomial::<_, u8>(&Q, a.get_vars().clone())
.to_number_field(&f);
let f2 = parse!("(1+a^2)+(a)*z+z^2")
.to_polynomial::<_, u8>(&Q, a.get_vars().clone())
.to_number_field(&f);
factors.sort_by(|a, b| a.0.internal_cmp(&b.0).then(a.1.cmp(&b.1)));
assert_eq!(factors, vec![(f1, 1), (f2, 1)])
}
}