From 768037586464e722b0c25d079ca6bd41711c40ff Mon Sep 17 00:00:00 2001 From: Ben Ruijl Date: Tue, 9 Jul 2024 20:43:59 +0200 Subject: [PATCH] Add factorization of polynomials over algebraic number fields --- src/domains/algebraic_number.rs | 86 ++++++++++++++++++++++++- src/poly/factor.rs | 108 +++++++++++++++++++++++++++++++- 2 files changed, 190 insertions(+), 4 deletions(-) diff --git a/src/domains/algebraic_number.rs b/src/domains/algebraic_number.rs index dd14ec2..880f5fd 100644 --- a/src/domains/algebraic_number.rs +++ b/src/domains/algebraic_number.rs @@ -5,7 +5,10 @@ use rand::Rng; use crate::{ coefficient::ConvertToRing, combinatorics::CombinationIterator, - poly::{factor::Factorize, gcd::PolynomialGCD, polynomial::MultivariatePolynomial, Variable}, + poly::{ + factor::Factorize, gcd::PolynomialGCD, polynomial::MultivariatePolynomial, Exponent, + Variable, + }, printer::PolynomialPrinter, }; @@ -541,6 +544,69 @@ impl> Field for AlgebraicExtension { } } +impl, E: Exponent> MultivariatePolynomial, E> { + /// Get the norm of a non-constant square-free polynomial `f` in the algebraic number field. + pub fn norm(&self) -> MultivariatePolynomial { + self.norm_impl().3 + } + + /// Get the norm of a non-constant square-free polynomial `f` in the algebraic number field. + /// Returns `(v, s, g, r)` where `v` is the shifted variable, `s` is the number of steps, + /// `g` is the shifted polynomial and `r` is the norm. + pub(crate) fn norm_impl( + &self, + ) -> ( + usize, + usize, + MultivariatePolynomial, + MultivariatePolynomial, + ) { + assert!(!self.is_constant()); + + let f = self.from_number_field(); + + let alpha = f + .get_vars() + .iter() + .position(|x| x == &self.field.poly.variables[0]) + .unwrap(); + + let mut poly = f.zero(); + let mut exp = vec![E::zero(); f.nvars()]; + for x in self.field.poly.into_iter() { + exp[alpha] = E::from_u32(x.exponents[0] as u32); + poly.append_monomial(x.coefficient.clone(), &exp); + } + + let poly_uni = poly.to_univariate(alpha); + + let mut s = 0; + loop { + let mut g_multi = f.clone(); + for v in 0..f.nvars() { + if v == alpha || f.degree(v) == E::zero() { + continue; + } + + let g_uni = g_multi.to_univariate(alpha); + + let r = g_uni.resultant_prs(&poly_uni); + + let d = r.derivative(v); + if r.gcd(&d).is_constant() { + return (v, s, g_multi, r); + } + + let alpha_poly = g_multi.variable(&self.get_vars_ref()[v]).unwrap() + - g_multi.variable(&self.field.poly.variables[0]).unwrap(); + + g_multi = g_multi.replace_with_poly(v, &alpha_poly); + s += 1; + } + } + } +} + #[cfg(test)] mod tests { use crate::atom::Atom; @@ -581,4 +647,22 @@ mod tests { } } } + + #[test] + fn norm() { + let a = Atom::parse("z^4+z^3+(2+a-a^2)z^2+(1+a^2-2a^3)z-2") + .unwrap() + .to_polynomial::<_, u8>(&Q, None); + let f = Atom::parse("a^4-3") + .unwrap() + .to_polynomial::<_, u16>(&Q, None); + let f = AlgebraicExtension::new(f); + let norm = a.to_number_field(&f).norm(); + + let res = Atom::parse("16-32*z-64*z^2-64*z^3-52*z^4-40*z^5-132*z^6-24*z^7-50*z^8+120*z^9+66*z^10+92*z^11+47*z^12+32*z^13+14*z^14+4*z^15+z^16") + .unwrap() + .to_polynomial::<_, u8>(&Q, a.variables.clone().into()); + + assert_eq!(norm, res); + } } diff --git a/src/poly/factor.rs b/src/poly/factor.rs index cfa90ab..1cd32c5 100644 --- a/src/poly/factor.rs +++ b/src/poly/factor.rs @@ -405,6 +405,80 @@ impl Factorize for MultivariatePolynomial Factorize + for MultivariatePolynomial, E, LexOrder> +{ + fn square_free_factorization(&self) -> Vec<(Self, usize)> { + if self.is_zero() { + return vec![]; + } + + let c = self.content(); + let stripped = self.clone().div_coeff(&c); + + let mut factors = vec![]; + + if !self.field.is_one(&c) { + factors.push((self.constant(c), 1)); + } + + 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() { + factors.push((self.one(), 1)) + } + + factors + } + + /// Perform Trager's algorithm for factorization. + fn factor(&self) -> Vec<(Self, usize)> { + let sf = self.square_free_factorization(); + + let mut full_factors = vec![]; + for (f, p) in &sf { + let (v, s, g, n) = f.norm_impl(); + + let factors = n.factor(); + + if factors.len() == 1 { + return vec![(f.clone(), 1)]; + } + + let mut g_f = g.to_number_field(&self.field); + + for (f, b) in factors { + debug!("Rational factor {}", f); + let alpha_poly = g.variable(&self.get_vars_ref()[v]).unwrap() + + g.variable(&self.field.poly().variables[0]).unwrap() + * &g.constant((s as u64).into()); + + let f = f.to_number_field(&self.field); + + let gcd = f.gcd(&g_f); + + g_f = g_f / &gcd; + + let g = MultivariatePolynomial::from_number_field(&gcd) + .replace_with_poly(v, &alpha_poly); + full_factors.push((g.to_number_field(&self.field), b * p)); + } + } + + full_factors + } + + fn is_irreducible(&self) -> bool { + // TODO: improve + self.factor().len() == 1 + } +} + impl< UField: FiniteFieldWorkspace, F: GaloisField> + PolynomialGCD, @@ -752,7 +826,7 @@ where for i in 0..2 * d { let r = self .field - .nth(rng.gen_range(0..self.field.characteristic().to_u64())); + .sample(&mut rng, (0, self.field.characteristic().to_u64() as i64)); if !F::is_zero(&r) { exp[var] = E::from_u32(i as u32); random_poly.append_monomial(r, &exp); @@ -770,8 +844,6 @@ where break g; } - // TODO: use Frobenius map and modular composition to prevent computing large exponent poly^(p^d) - let p: Integer = self.field.size().to_u64().into(); let b = if self.field.characteristic() == 2.into() { let max = self.field.get_extension_degree() as usize * d; @@ -785,6 +857,8 @@ where b } else { + // TODO: use Frobenius map and modular composition to prevent computing large exponent poly^(p^d) + let p = self.field.size(); random_poly .exp_mod_univariate(&(&p.pow(d as u64) - &1i64.into()) / &2i64.into(), &mut s) - self.one() @@ -3293,8 +3367,10 @@ mod test { use crate::{ atom::Atom, domains::{ + algebraic_number::AlgebraicExtension, finite_field::{Zp, Z2}, integer::Z, + rational::Q, }, poly::factor::Factorize, }; @@ -3532,4 +3608,30 @@ mod test { assert_eq!(a.factor().len(), 2) } + + #[test] + fn algebraic_extension() { + let a = Atom::parse("z^4+z^3+(2+a-a^2)z^2+(1+a^2-2a^3)z-2") + .unwrap() + .to_polynomial::<_, u8>(&Q, None); + let f = Atom::parse("a^4-3") + .unwrap() + .to_polynomial::<_, u16>(&Q, None); + let f = AlgebraicExtension::new(f); + + let mut factors = a.to_number_field(&f).factor(); + + let f1 = Atom::parse("(1-a^2)+(1-a)*z+z^2") + .unwrap() + .to_polynomial::<_, u8>(&Q, a.get_vars().clone().into()) + .to_number_field(&f); + let f2 = Atom::parse("(1+a^2)+(a)*z+z^2") + .unwrap() + .to_polynomial::<_, u8>(&Q, a.get_vars().clone().into()) + .to_number_field(&f); + + factors.sort_by(|a, b| a.partial_cmp(&b).unwrap()); + + assert_eq!(factors, vec![(f1, 1), (f2, 1)]) + } }