Add algebraic numbers and Galois fields to Python API

src/domains/algebraic_number.rs

@@ -1,4 +1,4 @@
-use std::{rc::Rc, sync::Arc};
+use std::sync::Arc;
 
 use rand::Rng;
 
@@ -25,7 +25,7 @@ use super::{
 // TODO: make special case for degree two and three and hardcode the multiplication table
 #[derive(Clone, PartialEq, Eq, PartialOrd, Hash)]
 pub struct AlgebraicExtension<R: Ring> {
-    poly: Rc<MultivariatePolynomial<R, u16>>, // TODO: convert to univariate polynomial
+    poly: Arc<MultivariatePolynomial<R, u16>>, // TODO: convert to univariate polynomial
 }
 
 impl<T: FiniteFieldWorkspace> GaloisField for AlgebraicExtension<FiniteField<T>>
@@ -62,7 +62,11 @@ where
         Self::Base: PolynomialGCD<u16>,
         <Self::Base as Ring>::Element: Copy,
     {
-        AlgebraicExtension::galois_field(self.poly.field.clone(), new_pow)
+        AlgebraicExtension::galois_field(
+            self.poly.field.clone(),
+            new_pow,
+            self.poly.variables[0].clone(),
+        )
     }
 
     fn upgrade_element(
@@ -155,12 +159,11 @@ where
 {
     /// Construct the Galois field GF(prime^exp).
     /// The irreducible polynomial is determined automatically.
-    pub fn galois_field(prime: FiniteField<UField>, exp: usize) -> Self {
+    pub fn galois_field(prime: FiniteField<UField>, exp: usize, var: Variable) -> Self {
         assert!(exp > 0);
 
         if exp == 1 {
-            let mut poly =
-                MultivariatePolynomial::new(&prime, None, Arc::new(vec![Variable::Temporary(0)]));
+            let mut poly = MultivariatePolynomial::new(&prime, None, Arc::new(vec![var]));
 
             poly.append_monomial(prime.one(), &[1]);
             return AlgebraicExtension::new(poly);
@@ -185,11 +188,7 @@ where
 
         let mut coeffs = vec![0; exp as usize + 1];
         coeffs[exp as usize] = 1;
-        let mut poly = MultivariatePolynomial::new(
-            &prime,
-            Some(coeffs.len()),
-            Arc::new(vec![Variable::Temporary(0)]),
-        );
+        let mut poly = MultivariatePolynomial::new(&prime, Some(coeffs.len()), Arc::new(vec![var]));
 
         // find the minimal polynomial
         let p = prime.get_prime().to_integer();
@@ -265,7 +264,7 @@ impl<R: Ring> AlgebraicExtension<R> {
     pub fn new(poly: MultivariatePolynomial<R, u16>) -> AlgebraicExtension<R> {
         if poly.nvars() == 1 {
             return AlgebraicExtension {
-                poly: Rc::new(poly),
+                poly: Arc::new(poly),
             };
         }
 
@@ -274,7 +273,7 @@ impl<R: Ring> AlgebraicExtension<R> {
         let uni = poly.to_univariate_from_univariate(v);
 
         AlgebraicExtension {
-            poly: Rc::new(uni.to_multivariate()),
+            poly: Arc::new(uni.to_multivariate()),
         }
     }
 
@@ -298,7 +297,7 @@ impl<R: Ring> AlgebraicExtension<R> {
         FiniteField<UField>: FiniteFieldCore<UField>,
     {
         AlgebraicExtension {
-            poly: Rc::new(
+            poly: Arc::new(
                 self.poly
                     .map_coeff(|c| c.to_finite_field(field), field.clone()),
             ),
@@ -663,6 +662,7 @@ mod tests {
     use crate::domains::algebraic_number::AlgebraicExtension;
     use crate::domains::finite_field::{PrimeIteratorU64, Zp, Z2};
     use crate::domains::rational::Q;
+    use crate::state::State;
 
     #[test]
     fn gcd_number_field() -> Result<(), String> {
@@ -688,12 +688,16 @@ mod tests {
     #[test]
     fn galois() {
         for j in 1..10 {
-            let _ = AlgebraicExtension::galois_field(Z2, j);
+            let _ = AlgebraicExtension::galois_field(Z2, j, State::get_symbol("v1").into());
         }
 
         for i in PrimeIteratorU64::new(2).take(20) {
             for j in 1..10 {
-                let _ = AlgebraicExtension::galois_field(Zp::new(i as u32), j);
+                let _ = AlgebraicExtension::galois_field(
+                    Zp::new(i as u32),
+                    j,
+                    State::get_symbol("v1").into(),
+                );
             }
         }
     }

src/domains/finite_field.rs

@@ -5,6 +5,7 @@ use std::ops::{Deref, Neg};
 
 use crate::domains::integer::Integer;
 use crate::poly::gcd::PolynomialGCD;
+use crate::poly::Variable::Temporary;
 use crate::printer::PrintOptions;
 
 use super::algebraic_number::AlgebraicExtension;
@@ -121,7 +122,7 @@ where
         Self::Base: PolynomialGCD<u16>,
         <Self::Base as Ring>::Element: Copy,
     {
-        AlgebraicExtension::galois_field(self.clone(), new_pow)
+        AlgebraicExtension::galois_field(self.clone(), new_pow, Temporary(0))
     }
 
     fn upgrade_element(

src/domains/rational.rs

@@ -12,7 +12,9 @@ use rug::{
 use crate::{poly::gcd::LARGE_U32_PRIMES, printer::PrintOptions, utils};
 
 use super::{
-    finite_field::{FiniteField, FiniteFieldCore, FiniteFieldWorkspace, ToFiniteField, Zp},
+    finite_field::{
+        FiniteField, FiniteFieldCore, FiniteFieldWorkspace, ToFiniteField, Two, Zp, Z2,
+    },
     integer::{Integer, Z},
     EuclideanDomain, Field, Ring,
 };
@@ -176,6 +178,26 @@ impl ToFiniteField<u32> for Rational {
     }
 }
 
+impl ToFiniteField<Two> for Rational {
+    fn to_finite_field(&self, field: &Z2) -> <Z2 as Ring>::Element {
+        match self {
+            &Rational::Natural(n, d) => {
+                let mut ff = n.rem_euclid(2) as u8;
+
+                if d != 1 {
+                    let df = n.rem_euclid(2) as u8;
+                    field.div_assign(&mut ff, &df);
+                }
+
+                ff
+            }
+            Rational::Large(r) => {
+                field.div(&(r.numer().mod_u(2) as u8), &(r.denom().mod_u(2) as u8))
+            }
+        }
+    }
+}
+
 impl Rational {
     pub fn new(mut num: i64, mut den: i64) -> Rational {
         if den == 0 {

src/poly.rs

@@ -1189,6 +1189,76 @@ impl<R: Ring, E: Exponent, O: MonomialOrder> MultivariatePolynomial<R, E, O> {
         out.as_view().normalize(workspace, &mut norm);
         std::mem::swap(norm.deref_mut(), out);
     }
+
+    pub fn to_expression_with_coeff_map<F: Fn(&R, &R::Element, &mut Atom)>(&self, f: F) -> Atom {
+        let mut out = Atom::default();
+        self.to_expression_with_coeff_map_into(f, &mut out);
+        out
+    }
+
+    pub fn to_expression_with_coeff_map_into<F: Fn(&R, &R::Element, &mut Atom)>(
+        &self,
+        f: F,
+        out: &mut Atom,
+    ) {
+        Workspace::get_local().with(|ws| self.to_expression_coeff_map_impl(ws, f, out));
+    }
+
+    pub(crate) fn to_expression_coeff_map_impl<F: Fn(&R, &R::Element, &mut Atom)>(
+        &self,
+        workspace: &Workspace,
+        f: F,
+        out: &mut Atom,
+    ) {
+        if self.is_zero() {
+            out.set_from_view(&workspace.new_num(0).as_view());
+            return;
+        }
+
+        let add = out.to_add();
+
+        let mut mul_h = workspace.new_atom();
+        let mut var_h = workspace.new_atom();
+        let mut num_h = workspace.new_atom();
+        let mut pow_h = workspace.new_atom();
+
+        let mut coeff = workspace.new_atom();
+        for monomial in self {
+            let mul = mul_h.to_mul();
+
+            for (var_id, &pow) in self.variables.iter().zip(monomial.exponents) {
+                if pow > E::zero() {
+                    match var_id {
+                        Variable::Symbol(v) => {
+                            var_h.to_var(*v);
+                        }
+                        Variable::Temporary(_) => {
+                            unreachable!("Temporary variables not supported");
+                        }
+                        Variable::Function(_, a) | Variable::Other(a) => {
+                            var_h.set_from_view(&a.as_view());
+                        }
+                    }
+
+                    if pow > E::one() {
+                        num_h.to_num((pow.to_u32() as i64).into());
+                        pow_h.to_pow(var_h.as_view(), num_h.as_view());
+                        mul.extend(pow_h.as_view());
+                    } else {
+                        mul.extend(var_h.as_view());
+                    }
+                }
+            }
+
+            f(&self.field, &monomial.coefficient, &mut coeff);
+            mul.extend(coeff.as_view());
+            add.extend(mul_h.as_view());
+        }
+
+        let mut norm = workspace.new_atom();
+        out.as_view().normalize(workspace, &mut norm);
+        std::mem::swap(norm.deref_mut(), out);
+    }
 }
 
 impl<R: Ring, E: Exponent> RationalPolynomial<R, E> {

src/poly/factor.rs

@@ -442,6 +442,11 @@ impl<E: Exponent> Factorize
 
         let mut full_factors = vec![];
         for (f, p) in &sf {
+            if f.is_constant() {
+                full_factors.push((f.clone(), *p));
+                continue;
+            }
+
             let (v, s, g, n) = f.norm_impl();
 
             let factors = n.factor();
@@ -811,12 +816,16 @@ where
         let mut exp = vec![E::zero(); self.nvars()];
 
         let mut try_counter = 0;
+        let characteristic = self.field.characteristic();
 
         let factor = loop {
             // generate a random non-constant polynomial
             random_poly.clear();
 
-            if d == 1 {
+            if d == 1
+                && (characteristic.is_zero()
+                    || &Integer::from(try_counter as i64) < &characteristic)
+            {
                 exp[var] = E::zero();
                 random_poly.append_monomial(self.field.nth(try_counter), &exp);
                 exp[var] = E::one();
@@ -826,7 +835,7 @@ where
                 for i in 0..2 * d {
                     let r = self
                         .field
-                        .sample(&mut rng, (0, self.field.characteristic().to_u64() as i64));
+                        .sample(&mut rng, (0, characteristic.to_i64().unwrap_or(i64::MAX)));
                     if !F::is_zero(&r) {
                         exp[var] = E::from_u32(i as u32);
                         random_poly.append_monomial(r, &exp);

src/poly/groebner.rs

@@ -3,7 +3,8 @@ use std::{cmp::Ordering, rc::Rc};
 use ahash::HashMap;
 
 use crate::domains::{
-    finite_field::{FiniteField, FiniteFieldCore, Mersenne64, Zp, Zp64},
+    algebraic_number::AlgebraicExtension,
+    finite_field::{FiniteField, FiniteFieldCore, Mersenne64, Zp, Zp64, Z2},
     rational::RationalField,
     Field, Ring,
 };
@@ -925,6 +926,10 @@ macro_rules! echelonize_impl {
 echelonize_impl!(Zp64);
 echelonize_impl!(FiniteField<Mersenne64>);
 echelonize_impl!(RationalField);
+echelonize_impl!(Z2);
+echelonize_impl!(AlgebraicExtension<Zp>);
+echelonize_impl!(AlgebraicExtension<Z2>);
+echelonize_impl!(AlgebraicExtension<RationalField>);
 
 #[cfg(test)]
 mod test {