alkahest_cas/poly/factor/
mod.rs

1//! V2-7 — Polynomial factorization over ℤ, 𝔽_p, and multivariate ℤ[𝑥₁,…].
2//!
3//! Univariate ℤ\[x\] uses FLINT `fmpz_poly_factor` (modular Berlekamp, Zassenhaus
4//! recombination, van Hoeij’s knapsack–LLL).  Multivariate ℤ[x₁,…] uses
5//! `fmpz_mpoly_factor` (Bernardin–Monagan EEZ pipeline).  Word-sized primes
6//! use `nmod_poly_factor` (Berlekamp / Cantor–Zassenhaus / Kaltofen–Shoup per
7//! FLINT’s internal choice).
8
9use super::error::FactorError;
10use super::multipoly::{multi_to_flint_pub, MultiPoly};
11use super::unipoly::UniPoly;
12use crate::flint::ffi::{self, FmpzMPolyFactorStruct, NmodPolyFactorStruct, NmodPolyStruct};
13use crate::flint::integer::FlintInteger;
14use crate::flint::mpoly::{FlintMPoly, FlintMPolyCtx};
15use crate::flint::FlintPoly;
16use crate::kernel::ExprId;
17
18#[cfg(not(flint3))]
19unsafe fn nmod_poly_factor_get_nth(
20    z: *mut NmodPolyStruct,
21    fac: *mut NmodPolyFactorStruct,
22    i: ffi::slong,
23) {
24    ffi::nmod_poly_factor_get_nmod_poly(z, fac, i);
25}
26
27#[cfg(flint3)]
28unsafe fn nmod_poly_factor_get_nth(
29    z: *mut NmodPolyStruct,
30    fac: *mut NmodPolyFactorStruct,
31    i: ffi::slong,
32) {
33    ffi::nmod_poly_factor_get_poly(z, fac as *const NmodPolyFactorStruct, i);
34}
35
36/// Factors of a non-zero `UniPoly`: `polynomial = unit · ∏ baseᵢ^expᵢ`.
37#[derive(Debug, Clone, PartialEq, Eq)]
38pub struct UniPolyFactorization {
39    pub unit: rug::Integer,
40    pub factors: Vec<(UniPoly, u32)>,
41}
42
43/// Factors of a non-zero multivariate integer polynomial.
44#[derive(Debug, Clone, PartialEq, Eq)]
45pub struct MultiPolyFactorization {
46    pub unit: rug::Integer,
47    pub factors: Vec<(MultiPoly, u32)>,
48}
49
50/// Factors of a univariate polynomial over ℤ/pℤ (coefficients ascending, reduced mod p).
51#[derive(Debug, Clone, PartialEq, Eq)]
52pub struct UniPolyFactorModP {
53    pub modulus: u64,
54    pub factors: Vec<(Vec<u64>, u32)>,
55}
56
57impl UniPolyFactorization {
58    /// Expand the factorization back to a single `UniPoly` (verification helper).
59    pub fn expand_with_var(&self, var: ExprId) -> UniPoly {
60        let mut acc = FlintPoly::from_rug_coefficients(std::slice::from_ref(&self.unit));
61        for (f, e) in &self.factors {
62            let powed = f.coeffs.pow(*e);
63            acc = &acc * &powed;
64        }
65        UniPoly { var, coeffs: acc }
66    }
67}
68
69impl MultiPolyFactorization {
70    /// Expand the factorization (verification helper).
71    pub fn expand_clone_vars(&self) -> MultiPoly {
72        let vars = self
73            .factors
74            .first()
75            .map(|(m, _)| m.vars.clone())
76            .unwrap_or_default();
77        let mut terms = std::collections::BTreeMap::new();
78        terms.insert(vec![], self.unit.clone());
79        let mut acc = MultiPoly { vars, terms };
80        for (f, e) in &self.factors {
81            let mut powered = MultiPoly::constant(f.vars.clone(), 1);
82            for _ in 0..*e {
83                powered = powered * f.clone();
84            }
85            acc = acc * powered;
86        }
87        acc
88    }
89}
90
91/// [`UniPoly::factor_z`](UniPoly::factor_z).
92pub fn factor_univariate_z(p: &UniPoly) -> Result<UniPolyFactorization, FactorError> {
93    p.factor_z()
94}
95
96/// [`MultiPoly::factor_z`](MultiPoly::factor_z).
97pub fn factor_multivariate_z(p: &MultiPoly) -> Result<MultiPolyFactorization, FactorError> {
98    p.factor_z()
99}
100
101impl UniPoly {
102    /// Factor over ℤ using FLINT (`fmpz_poly_factor`).
103    pub fn factor_z(&self) -> Result<UniPolyFactorization, FactorError> {
104        if self.is_zero() {
105            return Err(FactorError::ZeroPolynomial);
106        }
107        let (unit, facs) = self
108            .coeffs
109            .factor_over_z()
110            .map_err(|()| FactorError::FlintFailure)?;
111        let factors: Vec<_> = facs
112            .into_iter()
113            .map(|(c, e)| {
114                (
115                    UniPoly {
116                        var: self.var,
117                        coeffs: c,
118                    },
119                    e,
120                )
121            })
122            .collect();
123        Ok(UniPolyFactorization {
124            unit: unit.to_rug(),
125            factors,
126        })
127    }
128}
129
130impl MultiPoly {
131    /// Factor over ℤ[𝑥₁,…] using FLINT `fmpz_mpoly_factor`.
132    pub fn factor_z(&self) -> Result<MultiPolyFactorization, FactorError> {
133        if self.is_zero() {
134            return Err(FactorError::ZeroPolynomial);
135        }
136        let nvars = self.vars.len().max(1);
137        let ctx = FlintMPolyCtx::new(nvars);
138        let mut a = multi_to_flint_pub(self, &ctx);
139        unsafe {
140            let mut fac = std::mem::MaybeUninit::<FmpzMPolyFactorStruct>::uninit();
141            ffi::fmpz_mpoly_factor_init(fac.as_mut_ptr(), ctx.as_ptr());
142            let mut fac = fac.assume_init();
143            let ok = ffi::fmpz_mpoly_factor(&mut fac, a.as_ptr(), ctx.as_ptr());
144            if ok == 0 {
145                ffi::fmpz_mpoly_factor_clear(&mut fac, ctx.as_ptr());
146                a.clear_with_ctx(&ctx);
147                return Err(FactorError::FlintFailure);
148            }
149            if ffi::fmpz_cmp_ui(std::ptr::addr_of!(fac.constant_den), 1) != 0 {
150                ffi::fmpz_mpoly_factor_clear(&mut fac, ctx.as_ptr());
151                a.clear_with_ctx(&ctx);
152                return Err(FactorError::FlintFailure);
153            }
154            let mut unit = FlintInteger::new();
155            ffi::fmpz_mpoly_factor_get_constant_fmpz(unit.inner_mut_ptr(), &fac, ctx.as_ptr());
156            let n = ffi::fmpz_mpoly_factor_length(&fac, ctx.as_ptr());
157            let mut factors = Vec::with_capacity(n as usize);
158            for i in 0..n {
159                let mut base = FlintMPoly::new(&ctx);
160                ffi::fmpz_mpoly_factor_get_base(base.as_mut_ptr(), &fac, i, ctx.as_ptr());
161                let terms = base.terms(nvars, &ctx);
162                base.clear_with_ctx(&ctx);
163                let mp = MultiPoly {
164                    vars: self.vars.clone(),
165                    terms,
166                };
167                let exp = ffi::fmpz_mpoly_factor_get_exp_si(&mut fac, i, ctx.as_ptr()) as u32;
168                factors.push((mp, exp));
169            }
170            ffi::fmpz_mpoly_factor_clear(&mut fac, ctx.as_ptr());
171            a.clear_with_ctx(&ctx);
172            Ok(MultiPolyFactorization {
173                unit: unit.to_rug(),
174                factors,
175            })
176        }
177    }
178}
179
180/// Reduce coefficients mod `p` (must satisfy 2 ≤ p ≤ 2⁶³) and factor over 𝔽_p.
181pub fn factor_univariate_mod_p(
182    coeffs: &[i64],
183    modulus: u64,
184) -> Result<UniPolyFactorModP, FactorError> {
185    if modulus < 2 {
186        return Err(FactorError::InvalidModulus);
187    }
188    let p = modulus as i128;
189    let mut poly: NmodPolyStruct = unsafe { std::mem::zeroed() };
190    let factors = unsafe {
191        ffi::nmod_poly_init(&mut poly, modulus);
192        for (i, &c) in coeffs.iter().enumerate() {
193            let r = (c as i128 % p + p) % p;
194            ffi::nmod_poly_set_coeff_ui(&mut poly, i as ffi::slong, r as ffi::ulong);
195        }
196        let mut fac = std::mem::MaybeUninit::<NmodPolyFactorStruct>::uninit();
197        ffi::nmod_poly_factor_init(fac.as_mut_ptr());
198        let mut fac = fac.assume_init();
199        ffi::nmod_poly_factor(&mut fac, &poly);
200        let mut factors = Vec::with_capacity(fac.num as usize);
201        for i in 0..fac.num {
202            let mut z: NmodPolyStruct = std::mem::zeroed();
203            ffi::nmod_poly_init(&mut z, modulus);
204            nmod_poly_factor_get_nth(&mut z, &mut fac, i);
205            let deg = ffi::nmod_poly_degree(&z);
206            let mut vc = Vec::with_capacity(deg as usize + 1);
207            for j in 0..=deg {
208                vc.push(ffi::nmod_poly_get_coeff_ui(&z, j));
209            }
210            ffi::nmod_poly_clear(&mut z);
211            let exp = *fac.exp.add(i as usize) as u32;
212            factors.push((vc, exp));
213        }
214        ffi::nmod_poly_factor_clear(&mut fac);
215        ffi::nmod_poly_clear(&mut poly);
216        factors
217    };
218    Ok(UniPolyFactorModP { modulus, factors })
219}
220
221#[cfg(test)]
222mod tests {
223    use super::*;
224    use crate::kernel::{Domain, ExprPool};
225
226    #[test]
227    fn univariate_x_squared_minus_one() {
228        let pool = ExprPool::new();
229        let x = pool.symbol("x", Domain::Real);
230        let e = pool.add(vec![pool.pow(x, pool.integer(2_i32)), pool.integer(-1_i32)]);
231        let p = UniPoly::from_symbolic(e, x, &pool).unwrap();
232        let fac = p.factor_z().unwrap();
233        assert_eq!(fac.factors.len(), 2);
234        let prod = fac.expand_with_var(x);
235        assert_eq!(prod, p);
236    }
237
238    #[test]
239    fn swinnerton_dyer_irreducible_degree() {
240        let fp = FlintPoly::swinnerton_dyer(5);
241        assert_eq!(fp.degree(), 32);
242        let fac = fp.factor_over_z().unwrap();
243        assert_eq!(fac.1.len(), 1);
244        assert_eq!(fac.1[0].1, 1);
245    }
246
247    #[test]
248    fn cyclotomic_105_mod2_splits() {
249        let phi = FlintPoly::cyclotomic(105);
250        let coeffs: Vec<i64> = (0..phi.length())
251            .map(|i| {
252                (phi.get_coeff_flint(i).to_rug() % 2i32)
253                    .to_i64()
254                    .expect("coeff mod 2")
255            })
256            .collect();
257        let fac = factor_univariate_mod_p(&coeffs, 2).unwrap();
258        let deg_total: i64 = fac
259            .factors
260            .iter()
261            .map(|(f, e)| (f.len() as i64 - 1) * (*e as i64))
262            .sum();
263        assert_eq!(deg_total, phi.degree());
264        assert!(
265            fac.factors.len() >= 2,
266            "Φ_105 should have multiple factors over GF(2)"
267        );
268    }
269
270    #[test]
271    fn multivariate_product_recovered() {
272        let pool = ExprPool::new();
273        let x = pool.symbol("x", Domain::Real);
274        let y = pool.symbol("y", Domain::Real);
275        let vars = vec![x, y];
276        let f1 = MultiPoly::from_symbolic(
277            pool.add(vec![
278                pool.pow(x, pool.integer(2_i32)),
279                pool.pow(y, pool.integer(2_i32)),
280                pool.integer(-1_i32),
281            ]),
282            vars.clone(),
283            &pool,
284        )
285        .unwrap();
286        let x_minus_y = pool.add(vec![x, pool.mul(vec![pool.integer(-1i32), y])]);
287        let f2 = MultiPoly::from_symbolic(x_minus_y, vars.clone(), &pool).unwrap();
288        let product = f1.clone() * f2.clone();
289        let fac = product.factor_z().unwrap();
290        let expanded = fac.expand_clone_vars();
291        assert_eq!(expanded, product);
292    }
293}
alkahest_cas/poly/factor/mod.rs

alkahest_cas/poly/factor/
mod.rs
