Support merging of algebraic extensions

- Add missing power in resultant computation
benruijl · Jul 11, 2024 · a52797a · a52797a
1 parent 7680375
commit a52797a
Show file tree

Hide file tree

Showing 4 changed files with 128 additions and 23 deletions.
diff --git a/src/domains/algebraic_number.rs b/src/domains/algebraic_number.rs
@@ -260,9 +260,21 @@ where
 }
 
 impl<R: Ring> AlgebraicExtension<R> {
+    /// Create a new algebraic extension from a univariate polynomial.
+    /// The polynomial should be monic and irreducible.
     pub fn new(poly: MultivariatePolynomial<R, u16>) -> AlgebraicExtension<R> {
+        if poly.nvars() == 1 {
+            return AlgebraicExtension {
+                poly: Rc::new(poly),
+            };
+        }
+
+        assert_eq!((0..poly.nvars()).filter(|v| poly.degree(*v) > 0).count(), 1);
+        let v = (0..poly.nvars()).find(|v| poly.degree(*v) > 0).unwrap();
+        let uni = poly.to_univariate_from_univariate(v);
+
         AlgebraicExtension {
-            poly: Rc::new(poly),
+            poly: Rc::new(uni.to_multivariate()),
         }
     }
 
@@ -293,7 +305,7 @@ impl<R: Ring> AlgebraicExtension<R> {
         }
     }
 
-    pub fn to_element(&self, poly: MultivariatePolynomial<R, u16>) -> AlgebraicNumber<R> {
+    pub fn to_element(&self, poly: MultivariatePolynomial<R, u16>) -> <Self as Ring>::Element {
         assert!(poly.nvars() == 1);
 
         if poly.degree(0) >= self.poly.degree(0) {
@@ -544,6 +556,44 @@ impl<R: Field + PolynomialGCD<u16>> Field for AlgebraicExtension<R> {
     }
 }
 
+impl<R: Field + PolynomialGCD<u16>> AlgebraicExtension<R> {
+    /// Extend the current algebraic extension `R[a]` with `b`, whose minimal polynomial
+    /// is `R[a][b]` and form `R[b]`. Also return the new representation of `a` and `b`.
+    ///
+    /// `b`  must be irreducible over `R` and `R[a]`; this is not checked.
+    pub fn extend(
+        &self,
+        b: &MultivariatePolynomial<AlgebraicExtension<R>>,
+    ) -> (
+        AlgebraicExtension<R>,
+        <AlgebraicExtension<R> as Ring>::Element,
+        <AlgebraicExtension<R> as Ring>::Element,
+    )
+    where
+        AlgebraicExtension<R>: PolynomialGCD<u16>,
+        MultivariatePolynomial<R>: Factorize,
+        MultivariatePolynomial<AlgebraicExtension<R>>: Factorize,
+    {
+        assert_eq!(self, &b.field);
+
+        let (_, s, g, r) = b.norm_impl();
+        debug_assert!(r.is_irreducible());
+
+        let f: AlgebraicExtension<R> = AlgebraicExtension::new(r.clone());
+        let mut g2 = g.to_number_field(&f);
+        let mut h = self.poly.to_number_field(&f); // yields constant coeffs
+
+        g2.unify_variables(&mut h);
+        let g2 = g2.gcd(&h);
+
+        let a = f.neg(&f.div(&g2.get_constant(), &g2.lcoeff()));
+        let y = f.to_element(g2.field.poly.one().mul_exp(&[1]));
+        let b = f.sub(&y, &f.mul(&a, &f.nth(s as u64)));
+
+        (f, a, b)
+    }
+}
+
 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> {
@@ -582,12 +632,16 @@ impl<R: Field + PolynomialGCD<E>, E: Exponent> MultivariatePolynomial<AlgebraicE
 
         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;
                 }
 
+                // construct f(x-s*a)
+                let alpha_poly = f.variable(&self.get_vars_ref()[v]).unwrap()
+                    - f.variable(&self.field.poly.variables[0]).unwrap()
+                        * &f.constant(f.field.nth(s as u64));
+                let g_multi = f.clone().replace_with_poly(v, &alpha_poly);
                 let g_uni = g_multi.to_univariate(alpha);
 
                 let r = g_uni.resultant_prs(&poly_uni);
@@ -596,13 +650,9 @@ impl<R: Field + PolynomialGCD<E>, E: Exponent> MultivariatePolynomial<AlgebraicE
                 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;
             }
+
+            s += 1;
         }
     }
 }
@@ -665,4 +715,31 @@ mod tests {
 
         assert_eq!(norm, res);
     }
+
+    #[test]
+    fn extend() {
+        let a = Atom::parse("x^2-2").unwrap().to_polynomial(&Q, None);
+        let ae = AlgebraicExtension::new(a);
+
+        let b = Atom::parse("y^2-3")
+            .unwrap()
+            .to_polynomial(&Q, None)
+            .to_number_field(&ae);
+
+        let (c, rep1, rep2) = ae.extend(&b);
+
+        let rf = Atom::parse("1-10*y^2+y^4").unwrap().to_polynomial(&Q, None);
+
+        assert_eq!(c.poly.as_ref(), &rf);
+
+        let r1 = Atom::parse("-9/2y+1/2y^3")
+            .unwrap()
+            .to_polynomial::<_, u16>(&Q, None);
+        assert_eq!(rep1.poly, r1);
+
+        let r2 = Atom::parse("11/2*y-1/2*y^3")
+            .unwrap()
+            .to_polynomial::<_, u16>(&Q, None);
+        assert_eq!(rep2.poly, r2);
+    }
 }
diff --git a/src/poly/polynomial.rs b/src/poly/polynomial.rs
@@ -436,6 +436,15 @@ impl<F: Ring, E: Exponent, O: MonomialOrder> MultivariatePolynomial<F, E, O> {
         self.variables.as_ref()
     }
 
+    /// Renaname a variable.
+    pub fn rename_variable(&mut self, old: &Variable, new: &Variable) {
+        if let Some(pos) = self.variables.iter().position(|v| v == old) {
+            let mut new_vars = self.variables.as_ref().clone();
+            new_vars[pos] = new.clone();
+            self.variables = Arc::new(new_vars);
+        }
+    }
+
     /// Unify the variable maps of two polynomials, i.e.
     /// rewrite a polynomial in `x` and one in `y` to a
     /// two polynomial in `x` and `y`.

diff --git a/src/poly/resultant.rs b/src/poly/resultant.rs
@@ -9,6 +9,10 @@ impl<F: EuclideanDomain> UnivariatePolynomial<F> {
             return other.resultant_prs(self);
         }
 
+        if other.is_constant() {
+            return self.field.pow(&other.get_constant(), self.degree() as u64);
+        }
+
         let mut a = self.clone();
         let mut a_new = other.clone();
 
@@ -50,7 +54,7 @@ impl<F: EuclideanDomain> UnivariatePolynomial<F> {
             (a, a_new) = (a_new, r.div_coeff(&beta));
         }
 
-        a_new.coefficients.pop().unwrap()
+        self.field.pow(&a_new.coefficients.pop().unwrap(), deg)
     }
 }
 
@@ -122,18 +126,18 @@ mod test {
     #[test]
     fn resultant_prs_large() {
         let system = [
-        "-272853213601 + 114339252960*v2 - 4841413740*v2^2 + 296664007680*v4 - 25123011840*v2*v4 - 
-    32592015360*v4^2 - 4907531205*v5 + 6155208630*v5^2 - 3860046090*v5^3 + 1210353435*v5^4 - 
-    151807041*v5^5 + 312814245280*v6 - 97612876080*v2*v6 + 1518070410*v2^2*v6 - 
-    253265840640*v4*v6 + 7877554560*v2*v4*v6 + 10219530240*v4^2*v6 - 146051082720*v6^2 + 
-    29048482440*v2*v6^2 + 75369035520*v4*v6^2 + 35852138640*v6^3 - 3036140820*v2*v6^3 - 
+        "-272853213601 + 114339252960*v2 - 4841413740*v2^2 + 296664007680*v4 - 25123011840*v2*v4 -
+    32592015360*v4^2 - 4907531205*v5 + 6155208630*v5^2 - 3860046090*v5^3 + 1210353435*v5^4 -
+    151807041*v5^5 + 312814245280*v6 - 97612876080*v2*v6 + 1518070410*v2^2*v6 -
+    253265840640*v4*v6 + 7877554560*v2*v4*v6 + 10219530240*v4^2*v6 - 146051082720*v6^2 +
+    29048482440*v2*v6^2 + 75369035520*v4*v6^2 + 35852138640*v6^3 - 3036140820*v2*v6^3 -
     7877554560*v4*v6^3 - 4841413740*v6^4 + 303614082*v6^5",
-        "-121828703201 - 1128406464*v1 + 303614082*v1^2 + 24547177584*v2 - 303614082*v2^2 - 
-    2927757312*v3 + 1575510912*v1*v3 + 2043906048*v3^2 + 123022775808*v4 - 6600113280*v2*v4 - 
-    15080712192*v4^2 + 1480577211*v5 + 146055798*v5^2 - 1347744906*v5^3 + 816475707*v5^4 - 
-    151807041*v5^5 + 171636450272*v6 - 32717479104*v2*v6 + 303614082*v2^2*v6 - 
-    135541762560*v4*v6 + 3151021824*v2*v4*v6 + 6131718144*v4^2*v6 - 95005523376*v6^2 + 
-    13441077468*v2*v6^2 + 49947954048*v4*v6^2 + 26959925088*v6^3 - 1821684492*v2*v6^3 - 
+        "-121828703201 - 1128406464*v1 + 303614082*v1^2 + 24547177584*v2 - 303614082*v2^2 -
+    2927757312*v3 + 1575510912*v1*v3 + 2043906048*v3^2 + 123022775808*v4 - 6600113280*v2*v4 -
+    15080712192*v4^2 + 1480577211*v5 + 146055798*v5^2 - 1347744906*v5^3 + 816475707*v5^4 -
+    151807041*v5^5 + 171636450272*v6 - 32717479104*v2*v6 + 303614082*v2^2*v6 -
+    135541762560*v4*v6 + 3151021824*v2*v4*v6 + 6131718144*v4^2*v6 - 95005523376*v6^2 +
+    13441077468*v2*v6^2 + 49947954048*v4*v6^2 + 26959925088*v6^3 - 1821684492*v2*v6^3 -
     6302043648*v4*v6^3 - 4176745074*v6^4 + 303614082*v6^5",
     ];
 

diff --git a/src/poly/univariate.rs b/src/poly/univariate.rs
@@ -484,6 +484,21 @@ impl<F: Ring> UnivariatePolynomial<F> {
 
         res
     }
+
+    /// Convert from a univariate polynomial to a multivariate polynomial.
+    pub fn to_multivariate<E: Exponent>(self) -> MultivariatePolynomial<F, E> {
+        let mut res = MultivariatePolynomial::new(
+            &self.field,
+            self.degree().into(),
+            Arc::new(vec![self.variable.as_ref().clone()]),
+        );
+
+        for (p, c) in self.coefficients.into_iter().enumerate() {
+            res.append_monomial(c, &[E::from_u32(p as u32)]);
+        }
+
+        res
+    }
 }
 
 impl<F: Ring> PartialEq for UnivariatePolynomial<F> {
@@ -1016,8 +1031,8 @@ impl<F: Field> UnivariatePolynomial<F> {
 }
 
 impl<R: Ring, E: Exponent> UnivariatePolynomial<PolynomialRing<R, E>> {
-    // Convert from a univariate polynomial to a polynomial.
-    pub fn to_multivariate(self) -> MultivariatePolynomial<R, E> {
+    /// Convert a univariate polynomial of multivariate polynomials to a multivariate polynomial.
+    pub fn flatten(self) -> MultivariatePolynomial<R, E> {
         let Some(pos) = self
             .field
             .variables