Add factorization of polynomials over algebraic number fields

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<R: Field + PolynomialGCD<u16>> Field for AlgebraicExtension<R> {
     }
 }
 
+impl<R: Field + PolynomialGCD<E>, E: Exponent> MultivariatePolynomial<AlgebraicExtension<R>, E> {
+    /// Get the norm of a non-constant square-free polynomial `f` in the algebraic number field.
+    pub fn norm(&self) -> MultivariatePolynomial<R, E> {
+        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<R, E>,
+        MultivariatePolynomial<R, E>,
+    ) {
+        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);
+    }
 }

src/poly/factor.rs

@@ -405,6 +405,80 @@ impl<E: Exponent> Factorize for MultivariatePolynomial<RationalField, E, LexOrde
     }
 }
 
+impl<E: Exponent> Factorize
+    for MultivariatePolynomial<AlgebraicExtension<RationalField>, 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<Base = FiniteField<UField>> + PolynomialGCD<E>,
@@ -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)])
+    }
 }