From 3b5a2b022ae265dd8c5ef8a6b23bda532311439b Mon Sep 17 00:00:00 2001 From: Ben Ruijl Date: Mon, 8 Jul 2024 15:25:44 +0200 Subject: [PATCH] Abstract over prime fields and Galois field extensions - Factorization can now be performed over Galois fields - Automatically upgrade to larger Galois fields when field is too small for sampling during factorization - Add Z2 prime field with elements 0 and 1 - Make sure finite field to_integer does not truncate --- src/domains/algebraic_number.rs | 195 ++++++++++++--- src/domains/finite_field.rs | 408 ++++++++++++++++++++++++++++---- src/domains/integer.rs | 167 +++++++------ src/domains/rational.rs | 2 +- src/poly/factor.rs | 145 +++++++++--- src/poly/gcd.rs | 82 +------ src/poly/polynomial.rs | 2 +- 7 files changed, 742 insertions(+), 259 deletions(-) diff --git a/src/domains/algebraic_number.rs b/src/domains/algebraic_number.rs index a161d04..dd14ec2 100644 --- a/src/domains/algebraic_number.rs +++ b/src/domains/algebraic_number.rs @@ -3,14 +3,18 @@ use std::{rc::Rc, sync::Arc}; use rand::Rng; use crate::{ + coefficient::ConvertToRing, combinatorics::CombinationIterator, poly::{factor::Factorize, gcd::PolynomialGCD, polynomial::MultivariatePolynomial, Variable}, printer::PolynomialPrinter, }; use super::{ - finite_field::{FiniteField, FiniteFieldCore, FiniteFieldWorkspace, ToFiniteField}, + finite_field::{ + FiniteField, FiniteFieldCore, FiniteFieldWorkspace, GaloisField, ToFiniteField, + }, integer::Integer, + rational::Rational, EuclideanDomain, Field, Ring, }; @@ -18,38 +22,159 @@ 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 { - poly: Rc>, // TODO: convert to univariate polynomial + poly: Rc>, // TODO: convert to univariate polynomial +} + +impl GaloisField for AlgebraicExtension> +where + FiniteField: FiniteFieldCore + PolynomialGCD, +{ + type Base = FiniteField; + + fn get_extension_degree(&self) -> u64 { + self.poly.degree(0) as u64 + } + + fn to_integer(&self, a: &Self::Element) -> Integer { + let mut p = Integer::zero(); + for x in a.poly.into_iter() { + p = p + &(self.poly.field.to_integer(x.coefficient) + * &self.characteristic().pow(x.exponents[0] as u64)); + } + p + } + + fn to_symmetric_integer(&self, a: &Self::Element) -> Integer { + let r = self.to_integer(a); + let s = self.size(); + if &r * &2.into() > s { + &r - &s + } else { + r + } + } + + fn upgrade(&self, new_pow: usize) -> AlgebraicExtension + where + Self::Base: PolynomialGCD, + ::Element: Copy, + { + AlgebraicExtension::galois_field(self.poly.field.clone(), new_pow) + } + + fn upgrade_element( + &self, + e: &Self::Element, + larger_field: &AlgebraicExtension, + ) -> as Ring>::Element { + larger_field.to_element(e.poly.clone()) + } + + fn downgrade_element( + &self, + e: & as Ring>::Element, + ) -> Self::Element { + self.to_element(e.poly.clone()) + } +} + +impl ConvertToRing for AlgebraicExtension> +where + FiniteField: FiniteFieldCore + PolynomialGCD, + Integer: ToFiniteField, +{ + fn element_from_integer(&self, number: Integer) -> Self::Element { + let mut q = &number % &self.size(); + let mut pow = 0; + let mut poly = self.poly.zero(); + while !q.is_zero() { + let (qn, r) = q.quot_rem(&self.poly.field.size()); + poly.append_monomial(r.to_finite_field(&self.poly.field), &[pow]); + pow += 1; + q = qn; + } + + AlgebraicNumber { poly } + } + + fn element_from_coefficient(&self, number: crate::coefficient::Coefficient) -> Self::Element { + match number { + crate::coefficient::Coefficient::Rational(r) => { + let n = self.element_from_integer(r.numerator()); + let d = self.element_from_integer(r.denominator()); + self.div(&n, &d) + } + crate::coefficient::Coefficient::Float(_) => { + panic!("Cannot convert float coefficient to algebraic number") + } + crate::coefficient::Coefficient::FiniteField(_, _) => { + panic!("Cannot convert finite field coefficient to algebraic number") + } + crate::coefficient::Coefficient::RationalPolynomial(_) => { + panic!("Cannot convert rational polynomial coefficient to algebraic number") + } + } + } + + fn element_from_coefficient_view( + &self, + number: crate::coefficient::CoefficientView<'_>, + ) -> Self::Element { + match number { + crate::coefficient::CoefficientView::Natural(n, d) => { + let n = self.element_from_integer(n.into()); + let d = self.element_from_integer(d.into()); + self.div(&n, &d) + } + crate::coefficient::CoefficientView::Large(l) => { + let r = Rational::from_large(l.to_rat()); + let n = self.element_from_integer(r.numerator()); + let d = self.element_from_integer(r.denominator()); + self.div(&n, &d) + } + crate::coefficient::CoefficientView::Float(_) => { + panic!("Cannot convert float coefficient to algebraic number") + } + crate::coefficient::CoefficientView::FiniteField(_, _) => { + panic!("Cannot convert finite field coefficient to algebraic number") + } + crate::coefficient::CoefficientView::RationalPolynomial(_) => { + panic!("Cannot convert rational polynomial coefficient to algebraic number") + } + } + } } impl AlgebraicExtension> where - FiniteField: FiniteFieldCore + PolynomialGCD, + FiniteField: FiniteFieldCore + PolynomialGCD, + as Ring>::Element: Copy, { /// Construct the Galois field GF(prime^exp). /// The irreducible polynomial is determined automatically. - pub fn galois_field(prime: UField, exp: usize) -> Self { + pub fn galois_field(prime: FiniteField, exp: usize) -> Self { assert!(exp > 0); - let field = FiniteField::::new(prime.clone()); - if exp == 1 { let mut poly = - MultivariatePolynomial::new(&field, None, Arc::new(vec![Variable::Temporary(0)])); + MultivariatePolynomial::new(&prime, None, Arc::new(vec![Variable::Temporary(0)])); - poly.append_monomial(field.one(), &[1]); + poly.append_monomial(prime.one(), &[1]); return AlgebraicExtension::new(poly); } fn is_irreducible( coeffs: &[u64], - poly: &mut MultivariatePolynomial, u8>, + poly: &mut MultivariatePolynomial, u16>, ) -> bool where - FiniteField: FiniteFieldCore + PolynomialGCD, + FiniteField: FiniteFieldCore + PolynomialGCD, + as Ring>::Element: Copy, + AlgebraicExtension>: PolynomialGCD, { poly.clear(); for (i, c) in coeffs.iter().enumerate() { - poly.append_monomial(poly.field.nth(*c), &[i as u8]); + poly.append_monomial(poly.field.nth(*c), &[i as u16]); } poly.is_irreducible() @@ -58,13 +183,14 @@ where let mut coeffs = vec![0; exp as usize + 1]; coeffs[exp as usize] = 1; let mut poly = MultivariatePolynomial::new( - &field, + &prime, Some(coeffs.len()), Arc::new(vec![Variable::Temporary(0)]), ); // find the minimal polynomial - if prime.to_u64() == 2 { + let p = prime.get_prime().to_integer(); + if p == 2.into() { coeffs[0] = 1; // try all odd number of non-zero coefficients @@ -74,7 +200,7 @@ where let mut c = CombinationIterator::new(exp as usize - 1, g as usize); while let Some(comb) = c.next() { for i in 0..g as usize { - coeffs[comb[i]] = 1; + coeffs[comb[i] + 1] = 1; } if is_irreducible(&coeffs, &mut poly) { @@ -82,7 +208,7 @@ where } for i in 0..g as usize { - coeffs[comb[i]] = 0; + coeffs[comb[i] + 1] = 0; } } } @@ -90,8 +216,9 @@ where unreachable!("No irreducible polynomial found for GF({},{})", prime, exp); } + let sample_max = p.to_i64().unwrap_or(i64::MAX) as u64; if exp == 2 { - for k in 1..prime.to_u64() { + for k in 1..sample_max { coeffs[0] = k; if is_irreducible(&coeffs, &mut poly) { @@ -103,8 +230,8 @@ where } // try shape x^n+a*x+b for fast division - for k in 1..prime.to_u64() { - for k2 in 1..prime.to_u64() { + for k in 1..sample_max { + for k2 in 1..sample_max { coeffs[0] = k; coeffs[1] = k2; @@ -118,7 +245,7 @@ where let mut r = rand::thread_rng(); loop { for c in coeffs.iter_mut() { - *c = r.gen_range(0..prime.to_u64()); + *c = r.gen_range(0..sample_max); } coeffs[exp as usize] = 1; @@ -130,7 +257,7 @@ where } impl AlgebraicExtension { - pub fn new(poly: MultivariatePolynomial) -> AlgebraicExtension { + pub fn new(poly: MultivariatePolynomial) -> AlgebraicExtension { AlgebraicExtension { poly: Rc::new(poly), } @@ -143,7 +270,7 @@ impl AlgebraicExtension { } /// Get the minimal polynomial. - pub fn poly(&self) -> &MultivariatePolynomial { + pub fn poly(&self) -> &MultivariatePolynomial { &self.poly } @@ -163,7 +290,7 @@ impl AlgebraicExtension { } } - pub fn to_element(&self, poly: MultivariatePolynomial) -> AlgebraicNumber { + pub fn to_element(&self, poly: MultivariatePolynomial) -> AlgebraicNumber { assert!(poly.nvars() == 1); if poly.degree(0) >= self.poly.degree(0) { @@ -190,7 +317,7 @@ impl std::fmt::Display for AlgebraicExtension { #[derive(Clone, PartialEq, Eq, Hash)] pub struct AlgebraicNumber { - pub(crate) poly: MultivariatePolynomial, + pub(crate) poly: MultivariatePolynomial, } impl PartialOrd for AlgebraicNumber { @@ -227,7 +354,7 @@ impl AlgebraicNumber { } } - pub fn into_poly(self) -> MultivariatePolynomial { + pub fn into_poly(self) -> MultivariatePolynomial { self.poly } } @@ -322,10 +449,10 @@ impl Ring for AlgebraicExtension { } fn size(&self) -> Integer { - &self.poly.field.characteristic() * self.poly.degree(0) as i64 + self.poly.field.size().pow(self.poly.degree(0) as u64) } - /// Sample a monic polynomial. + /// Sample a polynomial. fn sample(&self, rng: &mut impl rand::RngCore, range: (i64, i64)) -> Self::Element { let coeffs: Vec<_> = (0..self.poly.degree(0)) .map(|_| self.poly.field.sample(rng, range)) @@ -334,7 +461,7 @@ impl Ring for AlgebraicExtension { let mut poly = self.poly.zero_with_capacity(coeffs.len()); let mut exp = vec![0]; for (i, c) in coeffs.into_iter().enumerate() { - exp[0] = i as u8; + exp[0] = i as u16; poly.append_monomial(c, &exp); } @@ -373,7 +500,7 @@ impl Ring for AlgebraicExtension { } } -impl> EuclideanDomain for AlgebraicExtension { +impl> EuclideanDomain for AlgebraicExtension { fn rem(&self, _a: &Self::Element, _b: &Self::Element) -> Self::Element { // TODO: due to the remainder requiring an inverse, we need to have R be a field // instead of a Euclidean domain. Relax this condition by doing a pseudo-division @@ -394,7 +521,7 @@ impl> EuclideanDomain for AlgebraicExtension { } } -impl> Field for AlgebraicExtension { +impl> Field for AlgebraicExtension { fn div(&self, a: &Self::Element, b: &Self::Element) -> Self::Element { self.mul(a, &self.inv(b)) } @@ -418,7 +545,7 @@ impl> Field for AlgebraicExtension { mod tests { use crate::atom::Atom; use crate::domains::algebraic_number::AlgebraicExtension; - use crate::domains::finite_field::{PrimeIteratorU64, Zp64}; + use crate::domains::finite_field::{PrimeIteratorU64, Zp, Z2}; use crate::domains::rational::Q; #[test] @@ -427,7 +554,7 @@ mod tests { let ring = AlgebraicExtension::new(ring); let a = Atom::parse("x^3-2x^2+(-2a^2+8a+2)x-a^2+11a-1")? - .to_polynomial::<_, u8>(&Q, None) + .to_polynomial::<_, u16>(&Q, None) .to_number_field(&ring); let b = Atom::parse("x^3-2x^2-x+1")? @@ -444,9 +571,13 @@ mod tests { #[test] fn galois() { + for j in 1..10 { + let _ = AlgebraicExtension::galois_field(Z2, j); + } + for i in PrimeIteratorU64::new(2).take(20) { for j in 1..10 { - let _ = AlgebraicExtension::::galois_field(i, j); + let _ = AlgebraicExtension::galois_field(Zp::new(i as u32), j); } } } diff --git a/src/domains/finite_field.rs b/src/domains/finite_field.rs index 1778e64..4238c7d 100644 --- a/src/domains/finite_field.rs +++ b/src/domains/finite_field.rs @@ -1,11 +1,13 @@ use rand::Rng; use std::fmt::{Display, Error, Formatter}; use std::hash::Hash; -use std::ops::Neg; +use std::ops::{Deref, Neg}; use crate::domains::integer::Integer; +use crate::poly::gcd::PolynomialGCD; use crate::printer::PrintOptions; +use super::algebraic_number::AlgebraicExtension; use super::integer::Z; use super::{EuclideanDomain, Field, Ring}; @@ -34,6 +36,110 @@ where ) -> as Ring>::Element; } +impl ToFiniteField for u32 { + fn to_finite_field(&self, field: &FiniteField) -> as Ring>::Element { + field.to_element(*self) + } +} + +impl ToFiniteField for u64 { + fn to_finite_field(&self, field: &FiniteField) -> as Ring>::Element { + field.to_element(*self) + } +} + +impl ToFiniteField for u32 { + fn to_finite_field(&self, field: &FiniteField) -> as Ring>::Element { + field.to_element(Two((*self % 2) as u8)) + } +} + +impl ToFiniteField for u64 { + fn to_finite_field(&self, field: &FiniteField) -> as Ring>::Element { + field.to_element(Two((*self % 2) as u8)) + } +} + +/// A Galois field `GF(p,n)` is a finite field with `p^n` elements. +/// It provides methods to upgrade and downgrade to Galois fields with the +/// same prime but with a different power. +pub trait GaloisField: Field { + type Base: Field; + + fn get_extension_degree(&self) -> u64; + // Convert a number from the finite field to standard form `[0,p)`. + fn to_integer(&self, a: &Self::Element) -> Integer; + /// Convert a number from the finite field to symmetric form `[-p/2,p/2]`. + fn to_symmetric_integer(&self, a: &Self::Element) -> Integer; + + /// Upgrade the field to `GF(p,new_pow)`. + fn upgrade(&self, new_pow: usize) -> AlgebraicExtension + where + Self::Base: PolynomialGCD, + ::Element: Copy; + + fn upgrade_element( + &self, + e: &Self::Element, + larger_field: &AlgebraicExtension, + ) -> as Ring>::Element; + + fn downgrade_element( + &self, + e: & as Ring>::Element, + ) -> Self::Element; +} + +impl GaloisField for FiniteField +where + FiniteField: Field + FiniteFieldCore, +{ + type Base = Self; + + fn get_extension_degree(&self) -> u64 { + 1 + } + + fn to_integer(&self, a: &Self::Element) -> Integer { + self.from_element(&a).to_integer() + } + + #[inline(always)] + fn to_symmetric_integer(&self, a: &Self::Element) -> Integer { + let i = self.from_element(a).to_integer(); + let p = self.get_prime().to_integer(); + + if &i * &2.into() > p { + &i - &p + } else { + i + } + } + + fn upgrade(&self, new_pow: usize) -> AlgebraicExtension> + where + Self::Base: PolynomialGCD, + ::Element: Copy, + { + AlgebraicExtension::galois_field(self.clone(), new_pow) + } + + fn upgrade_element( + &self, + e: &Self::Element, + larger_field: &AlgebraicExtension, + ) -> as Ring>::Element { + larger_field.constant(e.clone()) + } + + fn downgrade_element( + &self, + e: & as Ring>::Element, + ) -> Self::Element { + e.poly.get_constant() + } +} + /// A number in a finite field. #[derive(Debug, Copy, Clone, Hash, PartialEq, PartialOrd, Eq)] pub struct FiniteFieldElement(pub(crate) UField); @@ -41,6 +147,9 @@ pub struct FiniteFieldElement(pub(crate) UField); pub trait FiniteFieldWorkspace: Clone + Display + Eq + Hash { /// Convert to u64. fn to_u64(&self) -> u64; + fn to_integer(&self) -> Integer { + self.to_u64().into() + } } pub trait FiniteFieldCore: Field { @@ -50,19 +159,17 @@ pub trait FiniteFieldCore: Field { fn to_element(&self, a: UField) -> Self::Element; /// Convert a number from the finite field to standard form `[0,p)`. fn from_element(&self, a: &Self::Element) -> UField; - /// Convert a number from the finite field to symmetric form `[-p/2,p/2]`. - fn to_symmetric_integer(&self, a: &Self::Element) -> Integer; } /// The modular ring `Z / mZ`, where `m` can be any positive integer. In most cases, /// `m` will be a prime, and the domain will be a field. /// -/// `Zp` and `Zp64` use Montgomery modular arithmetic -/// to increase the performance of the multiplication operator. +/// [Zp] ([`FiniteField`]) and [Zp64] ([`FiniteField`]) use Montgomery modular arithmetic +/// to increase the performance of the multiplication operator. For the prime `2`, use [Z2] instead. /// -/// For `m` larger than `2^64`, use `FiniteField`. +/// For `m` larger than `2^64`, use [`FiniteField`]. /// -/// The special field `FiniteField` can be used to have even faster arithmetic +/// The special field [`FiniteField`] can be used to have even faster arithmetic /// for a field with Mersenne prime `2^61-1`. #[derive(Debug, Clone, PartialEq, Eq, Hash)] pub struct FiniteField { @@ -74,7 +181,9 @@ pub struct FiniteField { impl Zp { /// Create a new finite field. `n` must be a prime larger than 2. pub fn new(p: u32) -> Zp { - assert!(p % 2 != 0); + if p % 2 == 0 { + panic!("Prime 2 is not supported: use Z2 instead."); + } FiniteField { p, @@ -135,18 +244,6 @@ impl FiniteFieldCore for Zp { fn from_element(&self, a: &FiniteFieldElement) -> u32 { self.mul(a, &FiniteFieldElement(1)).0 } - - /// Convert a number from Montgomory form to symmetric form. - #[inline(always)] - fn to_symmetric_integer(&self, a: &FiniteFieldElement) -> Integer { - let i = self.from_element(a) as u64; - - if i * 2 > self.get_prime() as u64 { - &Integer::from(i) - &Integer::from(self.get_prime() as u64) - } else { - i.into() - } - } } impl Ring for Zp { @@ -379,7 +476,9 @@ impl FiniteFieldWorkspace for u64 { impl Zp64 { /// Create a new finite field. `n` must be a prime larger than 2. fn new(p: u64) -> Zp64 { - assert!(p % 2 != 0); + if p % 2 == 0 { + panic!("Prime 2 is not supported: use Z2 instead."); + } FiniteField { p, @@ -435,18 +534,6 @@ impl FiniteFieldCore for Zp64 { fn from_element(&self, a: &FiniteFieldElement) -> u64 { self.mul(a, &FiniteFieldElement(1)).0 } - - /// Convert a number from Montgomory form to symmetric form. - #[inline(always)] - fn to_symmetric_integer(&self, a: &FiniteFieldElement) -> Integer { - let i = self.from_element(a); - - if i > self.get_prime() / 2 { - &Integer::from(i) - &Integer::from(self.get_prime()) - } else { - i.into() - } - } } impl Display for FiniteField { @@ -669,6 +756,239 @@ impl Field for Zp64 { } } +/// The finite field with 0 and 1 as elements. +pub type Z2 = FiniteField; +/// The finite field with 0 and 1 as elements. +pub const Z2: FiniteField = Z2::new(); + +#[derive(Copy, Clone, Hash, Eq, PartialEq)] +pub struct Two(pub(crate) u8); + +impl Z2 { + /// Create a new finite field with prime 2. + pub const fn new() -> Z2 { + FiniteField { + p: Two(2), + m: Two(2), + one: FiniteFieldElement(Two(1)), + } + } + + // Get the prime 2. + pub fn get_prime() -> Two { + Two(2) + } +} + +impl Two { + pub const fn new() -> Two { + Two(2) + } +} + +impl Default for Two { + fn default() -> Self { + Two(2) + } +} + +impl Deref for Two { + type Target = u8; + + fn deref(&self) -> &Self::Target { + &self.0 + } +} + +impl std::fmt::Debug for Two { + fn fmt(&self, f: &mut Formatter<'_>) -> std::fmt::Result { + std::fmt::Debug::fmt(&self.0, f) + } +} + +impl Display for Two { + fn fmt(&self, f: &mut Formatter<'_>) -> std::fmt::Result { + self.0.fmt(f) + } +} + +impl FiniteFieldWorkspace for Two { + fn to_u64(&self) -> u64 { + self.0 as u64 + } +} + +impl FiniteFieldCore for FiniteField { + fn new(p: Two) -> Self { + FiniteField { + p, + m: p, + one: FiniteFieldElement(Two(1)), + } + } + + fn get_prime(&self) -> Two { + Two(2) + } + + fn to_element(&self, a: Two) -> Self::Element { + a.0 % 2 + } + + fn from_element(&self, a: &Self::Element) -> Two { + Two(*a) + } +} + +impl Ring for FiniteField { + type Element = u8; + + #[inline(always)] + fn add(&self, a: &Self::Element, b: &Self::Element) -> Self::Element { + a ^ b + } + + #[inline(always)] + fn sub(&self, a: &Self::Element, b: &Self::Element) -> Self::Element { + a ^ b + } + + #[inline(always)] + fn mul(&self, a: &Self::Element, b: &Self::Element) -> Self::Element { + *a * *b + } + + #[inline] + fn add_assign(&self, a: &mut Self::Element, b: &Self::Element) { + *a = self.add(a, b); + } + + #[inline] + fn sub_assign(&self, a: &mut Self::Element, b: &Self::Element) { + *a = self.sub(a, b); + } + + #[inline] + fn mul_assign(&self, a: &mut Self::Element, b: &Self::Element) { + *a = self.mul(a, b); + } + + fn add_mul_assign(&self, a: &mut Self::Element, b: &Self::Element, c: &Self::Element) { + self.add_assign(a, &self.mul(b, c)); + } + + fn sub_mul_assign(&self, a: &mut Self::Element, b: &Self::Element, c: &Self::Element) { + self.sub_assign(a, &self.mul(b, c)); + } + + /// Computes -x mod n. + #[inline] + fn neg(&self, a: &Self::Element) -> Self::Element { + *a + } + + #[inline] + fn zero(&self) -> Self::Element { + 0 + } + + /// Return the unit element in Montgomory form. + #[inline] + fn one(&self) -> Self::Element { + 1 + } + + #[inline] + fn nth(&self, n: u64) -> Self::Element { + (n % 2) as u8 + } + + /// Compute b^e % n. + #[inline] + fn pow(&self, b: &Self::Element, e: u64) -> Self::Element { + if e == 0 { + 1 + } else { + *b + } + } + + #[inline] + fn is_zero(a: &Self::Element) -> bool { + *a == 0 + } + + #[inline] + fn is_one(&self, a: &Self::Element) -> bool { + *a == 1 + } + + fn one_is_gcd_unit() -> bool { + true + } + + fn characteristic(&self) -> Integer { + 2.into() + } + + fn size(&self) -> Integer { + 2.into() + } + + fn sample(&self, rng: &mut impl rand::RngCore, _range: (i64, i64)) -> Self::Element { + rng.gen_range(0..2) + } + + fn fmt_display( + &self, + element: &Self::Element, + opts: &PrintOptions, + _in_product: bool, + f: &mut Formatter<'_>, + ) -> Result<(), Error> { + if opts.symmetric_representation_for_finite_field { + self.to_symmetric_integer(element).fmt(f) + } else { + self.from_element(element).fmt(f) + } + } +} + +impl EuclideanDomain for FiniteField { + #[inline] + fn rem(&self, _: &Self::Element, _: &Self::Element) -> Self::Element { + 0 + } + + #[inline] + fn quot_rem(&self, a: &Self::Element, b: &Self::Element) -> (Self::Element, Self::Element) { + (self.mul(a, &self.inv(b)), 0) + } + + #[inline] + fn gcd(&self, _: &Self::Element, _: &Self::Element) -> Self::Element { + 1 + } +} + +impl Field for FiniteField { + #[inline] + fn div(&self, a: &Self::Element, b: &Self::Element) -> Self::Element { + self.mul(a, &self.inv(b)) + } + + #[inline] + fn div_assign(&self, a: &mut Self::Element, b: &Self::Element) { + *a = self.mul(a, &self.inv(b)); + } + + /// Computes x^-1 mod n. + fn inv(&self, a: &Self::Element) -> Self::Element { + assert!(*a != 0, "0 is not invertible"); + 1 + } +} + /// The 64-bit Mersenne prime 2^61 -1. /// /// Can be used for faster finite field arithmetic @@ -733,14 +1053,6 @@ impl FiniteFieldCore for FiniteField { fn from_element(&self, a: &Self::Element) -> Mersenne64 { Mersenne64(*a) } - - fn to_symmetric_integer(&self, a: &Self::Element) -> Integer { - if *a * 2 > Mersenne64::PRIME { - &Integer::from(*a) - &Integer::from(Mersenne64::PRIME) - } else { - (*a).into() - } - } } impl Ring for FiniteField { @@ -957,6 +1269,10 @@ impl FiniteFieldWorkspace for Integer { panic!("Modulus {} is too large to be converted to u64", self) } } + + fn to_integer(&self) -> Integer { + self.clone() + } } /// A finite field with a large prime modulus. @@ -980,11 +1296,11 @@ impl FiniteFieldCore for FiniteField { } fn from_element(&self, a: &Integer) -> Integer { - a.clone() - } - - fn to_symmetric_integer(&self, a: &Integer) -> Integer { - a.clone() + if a.is_negative() { + a.clone() + &self.p + } else { + a.clone() + } } } diff --git a/src/domains/integer.rs b/src/domains/integer.rs index bc2fb82..1ce27e4 100644 --- a/src/domains/integer.rs +++ b/src/domains/integer.rs @@ -16,7 +16,8 @@ use crate::{printer::PrintOptions, utils}; use super::{ finite_field::{ - FiniteField, FiniteFieldCore, FiniteFieldWorkspace, Mersenne64, ToFiniteField, Zp, Zp64, + FiniteField, FiniteFieldCore, FiniteFieldWorkspace, Mersenne64, ToFiniteField, Two, Zp, + Zp64, Z2, }, rational::Rational, EuclideanDomain, Ring, @@ -168,6 +169,16 @@ impl ToFiniteField for Integer { } } +impl ToFiniteField for Integer { + fn to_finite_field(&self, field: &Z2) -> ::Element { + match self { + &Integer::Natural(n) => field.to_element(Two(n.rem_euclid(2 as i64) as u8)), + &Integer::Double(n) => field.to_element(Two(n.rem_euclid(2 as i128) as u8)), + Integer::Large(r) => field.to_element(Two(r.mod_u(2) as u8)), + } + } +} + impl ToFiniteField for Integer { fn to_finite_field( &self, @@ -473,6 +484,84 @@ impl Integer { } } + pub fn quot_rem(&self, b: &Integer) -> (Integer, Integer) { + if b.is_zero() { + panic!("Cannot divide by zero"); + } + + match (self, b) { + (Integer::Natural(aa), Integer::Natural(bb)) => { + if let Some(q) = aa.checked_div_euclid(*bb) { + (Integer::Natural(q), self - &(b * &Integer::Natural(q))) + } else { + (Integer::Double(-(i64::MIN as i128)), Integer::zero()) + } + } + (Integer::Natural(a), Integer::Double(b)) => { + // we always have |a| <= |b| + if *a < 0 { + if *b > 0 { + (Integer::Natural(-1), Integer::from_double(*a as i128 + *b)) + } else { + (Integer::Natural(1), Integer::from_double(*a as i128 - *b)) + } + } else { + (Integer::zero(), Integer::Natural(*a)) + } + } + (Integer::Double(aa), Integer::Natural(bb)) => { + if let Some(q) = aa.checked_div_euclid(*bb as i128) { + let q = Integer::from_double(q); + (q.clone(), self - &(b * &q)) + } else { + ( + Integer::Large(MultiPrecisionInteger::from(i128::MIN).neg()), + Integer::zero(), + ) + } + } + (Integer::Double(aa), Integer::Double(bb)) => { + let q = Integer::from_double(aa.div_euclid(*bb)); // b != -1 + (q.clone(), self - &(b * &q)) + } + (Integer::Natural(a), Integer::Large(b)) => { + if *a < 0 { + if *b > 0 { + (Integer::Natural(-1), Integer::from_large((a + b).into())) + } else { + (Integer::Natural(1), Integer::from_large((a - b).into())) + } + } else { + (Integer::zero(), Integer::Natural(*a)) + } + } + (Integer::Large(a), Integer::Natural(b)) => { + let r = a.clone().div_rem_euc(MultiPrecisionInteger::from(*b)); + (Integer::from_large(r.0), Integer::from_large(r.1)) + } + (Integer::Large(a), Integer::Large(b)) => { + let r = a.clone().div_rem_euc(b.clone()); + (Integer::from_large(r.0), Integer::from_large(r.1)) + } + + (Integer::Double(a), Integer::Large(b)) => { + if *a < 0 { + if *b > 0 { + (Integer::Natural(-1), Integer::from_large((a + b).into())) + } else { + (Integer::Natural(1), Integer::from_large((a - b).into())) + } + } else { + (Integer::zero(), Integer::Double(*a)) + } + } + (Integer::Large(a), Integer::Double(b)) => { + let r = a.clone().div_rem_euc(MultiPrecisionInteger::from(*b)); + (Integer::from_large(r.0), Integer::from_large(r.1)) + } + } + } + pub fn gcd(&self, b: &Integer) -> Integer { match (self, b) { (Integer::Natural(n1), Integer::Natural(n2)) => { @@ -831,81 +920,7 @@ impl EuclideanDomain for IntegerRing { } fn quot_rem(&self, a: &Self::Element, b: &Self::Element) -> (Self::Element, Self::Element) { - if b.is_zero() { - panic!("Cannot divide by zero"); - } - - match (a, b) { - (Integer::Natural(aa), Integer::Natural(bb)) => { - if let Some(q) = aa.checked_div_euclid(*bb) { - (Integer::Natural(q), a - &(b * &Integer::Natural(q))) - } else { - (Integer::Double(-(i64::MIN as i128)), Integer::zero()) - } - } - (Integer::Natural(a), Integer::Double(b)) => { - // we always have |a| <= |b| - if *a < 0 { - if *b > 0 { - (Integer::Natural(-1), Integer::from_double(*a as i128 + *b)) - } else { - (Integer::Natural(1), Integer::from_double(*a as i128 - *b)) - } - } else { - (Integer::zero(), Integer::Natural(*a)) - } - } - (Integer::Double(aa), Integer::Natural(bb)) => { - if let Some(q) = aa.checked_div_euclid(*bb as i128) { - let q = Integer::from_double(q); - (q.clone(), a - &(b * &q)) - } else { - ( - Integer::Large(MultiPrecisionInteger::from(i128::MIN).neg()), - Integer::zero(), - ) - } - } - (Integer::Double(aa), Integer::Double(bb)) => { - let q = Integer::from_double(aa.div_euclid(*bb)); // b != -1 - (q.clone(), a - &(b * &q)) - } - (Integer::Natural(a), Integer::Large(b)) => { - if *a < 0 { - if *b > 0 { - (Integer::Natural(-1), Integer::from_large((a + b).into())) - } else { - (Integer::Natural(1), Integer::from_large((a - b).into())) - } - } else { - (Integer::zero(), Integer::Natural(*a)) - } - } - (Integer::Large(a), Integer::Natural(b)) => { - let r = a.clone().div_rem_euc(MultiPrecisionInteger::from(*b)); - (Integer::from_large(r.0), Integer::from_large(r.1)) - } - (Integer::Large(a), Integer::Large(b)) => { - let r = a.clone().div_rem_euc(b.clone()); - (Integer::from_large(r.0), Integer::from_large(r.1)) - } - - (Integer::Double(a), Integer::Large(b)) => { - if *a < 0 { - if *b > 0 { - (Integer::Natural(-1), Integer::from_large((a + b).into())) - } else { - (Integer::Natural(1), Integer::from_large((a - b).into())) - } - } else { - (Integer::zero(), Integer::Double(*a)) - } - } - (Integer::Large(a), Integer::Double(b)) => { - let r = a.clone().div_rem_euc(MultiPrecisionInteger::from(*b)); - (Integer::from_large(r.0), Integer::from_large(r.1)) - } - } + a.quot_rem(b) } fn gcd(&self, a: &Self::Element, b: &Self::Element) -> Self::Element { diff --git a/src/domains/rational.rs b/src/domains/rational.rs index 5e82e75..9e260ca 100644 --- a/src/domains/rational.rs +++ b/src/domains/rational.rs @@ -553,7 +553,7 @@ impl Rational { let p = LARGE_U32_PRIMES[prime_start]; // TODO: support u64 prime_start += 1; - let field = FiniteField::new(p); + let field = FiniteField::::new(p); prime_sample_point.clear(); for x in sample { prime_sample_point.push(x.to_finite_field(&field)); diff --git a/src/poly/factor.rs b/src/poly/factor.rs index 6e028d7..cfa90ab 100644 --- a/src/poly/factor.rs +++ b/src/poly/factor.rs @@ -7,9 +7,10 @@ use tracing::debug; use crate::{ combinatorics::CombinationIterator, domains::{ + algebraic_number::AlgebraicExtension, finite_field::{ - FiniteField, FiniteFieldCore, FiniteFieldWorkspace, PrimeIteratorU64, ToFiniteField, - Zp, Zp64, + FiniteField, FiniteFieldCore, FiniteFieldWorkspace, GaloisField, PrimeIteratorU64, + ToFiniteField, Zp, Zp64, }, integer::{Integer, IntegerRing, Z}, rational::{RationalField, Q}, @@ -404,10 +405,15 @@ impl Factorize for MultivariatePolynomial Factorize - for MultivariatePolynomial, E, LexOrder> +impl< + UField: FiniteFieldWorkspace, + F: GaloisField> + PolynomialGCD, + E: Exponent, + > Factorize for MultivariatePolynomial where - FiniteField: Field + PolynomialGCD + FiniteFieldCore, + FiniteField: Field + FiniteFieldCore + PolynomialGCD, + as Ring>::Element: Copy, + AlgebraicExtension<::Base>: PolynomialGCD, { fn square_free_factorization(&self) -> Vec<(Self, usize)> { let c = self.lcoeff(); @@ -554,10 +560,15 @@ where } } -impl - MultivariatePolynomial, E, LexOrder> +impl< + UField: FiniteFieldWorkspace, + F: GaloisField> + PolynomialGCD, + E: Exponent, + > MultivariatePolynomial where - FiniteField: Field + PolynomialGCD + FiniteFieldCore, + FiniteField: Field + FiniteFieldCore + PolynomialGCD, + as Ring>::Element: Copy, + AlgebraicExtension<::Base>: PolynomialGCD, { /// Bernardin's algorithm for square free factorization. fn square_free_factorization_bernardin(&self) -> Vec<(Self, usize)> { @@ -590,7 +601,7 @@ where // take the pth root // the coefficients remain unchanged, since x^1/p = x // since the derivative in every var is 0, all powers are divisible by p - let p = self.field.get_prime().to_u64() as usize; + let p = self.field.characteristic().to_u64() as usize; let mut b = f.clone(); for es in b.exponents_iter_mut() { for e in es { @@ -648,7 +659,7 @@ where let mut factors = vec![]; let mut i = 1; - while !w.is_constant() && i < self.field.get_prime().to_u64() as usize { + while !w.is_constant() && i < self.field.characteristic().to_u64() as usize { let z = v - w.derivative(var); let g = w.gcd(&z); w = w / &g; @@ -666,7 +677,6 @@ where /// Perform distinct degree factorization on a monic, univariate and square-free polynomial. pub fn distinct_degree_factorization(&self) -> Vec<(usize, Self)> { - assert!(self.field.get_prime().to_u64() != 2); let Some(var) = self.last_exponents().iter().position(|x| *x > E::zero()) else { return vec![(0, self.clone())]; // constant polynomial }; @@ -682,7 +692,7 @@ where while !f.is_one() { i += 1; - h = h.exp_mod_univariate(self.field.get_prime().to_u64().into(), &mut f); + h = h.exp_mod_univariate(self.field.size(), &mut f); let mut g = f.gcd(&(&h - &x)); @@ -706,7 +716,6 @@ where /// Perform Cantor-Zassenhaus's probabilistic algorithm for /// finding irreducible factors of degree `d`. pub fn equal_degree_factorization(&self, d: usize) -> Vec { - assert!(self.field.get_prime().to_u64() != 2); let mut s = self.clone().make_monic(); let Some(var) = self.last_exponents().iter().position(|x| *x > E::zero()) else { @@ -743,8 +752,8 @@ where for i in 0..2 * d { let r = self .field - .nth(rng.gen_range(0..self.field.get_prime().to_u64())); - if !FiniteField::::is_zero(&r) { + .nth(rng.gen_range(0..self.field.characteristic().to_u64())); + if !F::is_zero(&r) { exp[var] = E::from_u32(i as u32); random_poly.append_monomial(r, &exp); } @@ -762,10 +771,24 @@ where } // TODO: use Frobenius map and modular composition to prevent computing large exponent poly^(p^d) - let p: Integer = self.field.get_prime().to_u64().into(); - let b = random_poly - .exp_mod_univariate(&(&p.pow(d as u64) - &1i64.into()) / &2i64.into(), &mut s) - - self.one(); + 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; + + let mut b = random_poly.clone(); + let mut vcur = b.clone(); + + for _ in 1..max { + vcur = (&vcur * &vcur).rem(&s); + b = b + vcur.clone(); + } + + b + } else { + random_poly + .exp_mod_univariate(&(&p.pow(d as u64) - &1i64.into()) / &2i64.into(), &mut s) + - self.one() + }; let g = b.gcd(&s); @@ -910,21 +933,43 @@ where let mut uni_f = self.replace(interpolation_var, &sample_point); let mut i = 0; + let mut rng = thread_rng(); loop { + if Integer::from(i) == self.field.size() { + let field = self + .field + .upgrade(self.field.get_extension_degree().to_u64() as usize + 1); + + debug!( + "Upgrading to Galois field with exponent {}", + field.get_extension_degree() + ); + + let s_l = self.map_coeff(|c| self.field.upgrade_element(c, &field), field.clone()); + + let facs = s_l.bivariate_factorization(main_var, interpolation_var); + + return facs + .into_iter() + .map(|f| f.map_coeff(|c| self.field.downgrade_element(c), self.field.clone())) + .collect(); + } + if self.degree(main_var) == uni_f.degree(main_var) && uni_f.gcd(&uni_f.derivative(main_var)).is_constant() { break; } - sample_point = self.field.nth(i); + // TODO: sample simple points first + sample_point = self.field.sample(&mut rng, (0, i)); uni_f = self.replace(interpolation_var, &sample_point); i += 1; } let mut d = self.degree(interpolation_var).to_u32(); - let shifted_poly = if !FiniteField::::is_zero(&sample_point) { + let shifted_poly = if !F::is_zero(&sample_point) { self.shift_var_cached(interpolation_var, &sample_point) } else { self.clone() @@ -993,7 +1038,7 @@ where rec_factors.push(rest); - if !FiniteField::::is_zero(&sample_point) { + if !F::is_zero(&sample_point) { for x in &mut rec_factors { // shift the polynomial to y - sample *x = x.shift_var_cached(interpolation_var, &self.field.neg(&sample_point)); @@ -1022,8 +1067,8 @@ where fn canonical_sort( biv_polys: &[Self], replace_var: usize, - sample_points: &[(usize, as Ring>::Element)], - ) -> Vec<(Self, as Ring>::Element, Self)> { + sample_points: &[(usize, ::Element)], + ) -> Vec<(Self, ::Element, Self)> { let mut univariate_factors = biv_polys .iter() .map(|f| { @@ -1051,7 +1096,7 @@ where fn lcoeff_precomputation( &self, bivariate_factors: &[Self], - sample_points: &[(usize, as Ring>::Element)], + sample_points: &[(usize, ::Element)], order: &[usize], ) -> Result<(Vec, Vec), usize> { let lcoeff = self.univariate_lcoeff(order[0]); @@ -1218,7 +1263,7 @@ where fn multivariate_hensel_lift_with_auto_lcoeff_fixing( &self, factors: &[Self], - sample_points: &[(usize, as Ring>::Element)], + sample_points: &[(usize, ::Element)], order: &[usize], ) -> Vec { let lcoeff = self.univariate_lcoeff(order[0]); @@ -1282,7 +1327,7 @@ where fn univariate_diophantine_field( factors: &[Self], order: &[usize], - sample_points: &[(usize, as Ring>::Element)], + sample_points: &[(usize, ::Element)], ) -> (Vec, Vec) { // produce univariate factors and univariate delta let mut univariate_factors = factors.to_vec(); @@ -1351,7 +1396,35 @@ where let uni_lcoeff = self.univariate_lcoeff(order[0]); let mut content_fail_count = 0; + let mut sample_fail = Integer::zero(); 'new_sample: loop { + sample_fail = sample_fail + &1.into(); + + if &sample_fail * &2.into() > self.field.size() { + // the field is too small, upgrade + let field = self + .field + .upgrade(self.field.get_extension_degree().to_u64() as usize + 1); + + debug!( + "Upgrading to Galois field with exponent {}", + field.get_extension_degree() + ); + + let s_l = self.map_coeff(|c| self.field.upgrade_element(c, &field), field.clone()); + + let facs = s_l.multivariate_factorization( + order, + coefficient_upper_bound, + max_bivariate_factors, + ); + + return facs + .into_iter() + .map(|f| f.map_coeff(|c| self.field.downgrade_element(c), self.field.clone())) + .collect(); + } + for s in &mut sample_points { s.1 = self.field.nth(rng.gen_range(0..=coefficient_upper_bound)); } @@ -1413,7 +1486,7 @@ where } for (v, s) in &sample_points { - debug!("Sample point {}={}", v, self.field.from_element(s)); + debug!("Sample point {}={}", v, self.field.printer(s)); } let bivariate_factors = biv_f.bivariate_factorization(order[0], order[1]); @@ -3219,7 +3292,10 @@ impl MultivariatePolynomial, E, LexOrder> { mod test { use crate::{ atom::Atom, - domains::{finite_field::Zp, integer::Z}, + domains::{ + finite_field::{Zp, Z2}, + integer::Z, + }, poly::factor::Factorize, }; @@ -3445,4 +3521,15 @@ mod test { r.sort_by(|a, b| a.partial_cmp(&b).unwrap()); assert_eq!(r, res); } + + #[test] + fn galois_upgrade() { + let a = Atom::parse( + "x^7(y^5+y^4+y^3+y^2)+x^5(y^3+y)+x^4(y^4+y)+x^3(y^2+y)+x^2y+x*y^2+x*y+x+y+1", + ) + .unwrap() + .to_polynomial::<_, u8>(&Z2, None); + + assert_eq!(a.factor().len(), 2) + } } diff --git a/src/poly/gcd.rs b/src/poly/gcd.rs index ee5d2cc..d566cd6 100755 --- a/src/poly/gcd.rs +++ b/src/poly/gcd.rs @@ -9,7 +9,7 @@ use tracing::{debug, instrument}; use crate::domains::algebraic_number::AlgebraicExtension; use crate::domains::finite_field::{ - FiniteField, FiniteFieldCore, FiniteFieldWorkspace, ToFiniteField, Zp, + FiniteField, FiniteFieldCore, FiniteFieldWorkspace, GaloisField, ToFiniteField, Zp, }; use crate::domains::integer::{FromFiniteField, Integer, IntegerRing, SMALL_PRIMES, Z}; use crate::domains::rational::{Rational, RationalField, Q}; @@ -2610,7 +2610,8 @@ impl PolynomialGCD for RationalField { } } -impl PolynomialGCD for FiniteField +impl>, E: Exponent> + PolynomialGCD for F where FiniteField: FiniteFieldCore, as Ring>::Element: Copy, @@ -2640,26 +2641,11 @@ where // upgrade to a Galois field that is large enough // TODO: start at a better bound? // TODO: run with Zp[var]/m_i instead and use CRT - for i in 2..100 { - debug!("Upgrading to Galois field: {}^{}", a.field.get_prime(), i); - - let field = AlgebraicExtension::galois_field(a.field.get_prime(), i); - - let ag = a.map_coeff(|c| field.constant(c.clone()), field.clone()); - let bg = b.map_coeff(|c| field.constant(c.clone()), field.clone()); - - if let Some(g) = MultivariatePolynomial::gcd_shape_modular( - &ag, - &bg, - vars, - bounds, - tight_bounds, - ) { - return g.map_coeff(|c| c.poly.get_constant(), a.field.clone()); - } - } - - panic!("Failed to find a suitable Galois field for gcd"); + let field = a.field.upgrade(a.field.get_extension_degree() as usize + 1); + let ag = a.map_coeff(|c| a.field.upgrade_element(c, &field), field.clone()); + let bg = b.map_coeff(|c| a.field.upgrade_element(c, &field), field.clone()); + let g = PolynomialGCD::gcd(&ag, &bg, vars, bounds, tight_bounds); + g.map_coeff(|c| a.field.downgrade_element(c), a.field.clone()) } } } @@ -3081,55 +3067,3 @@ impl PolynomialGCD for AlgebraicExtension { } } } - -impl PolynomialGCD - for AlgebraicExtension> -where - FiniteField: FiniteFieldCore, - as Ring>::Element: Copy, -{ - fn heuristic_gcd( - _a: &MultivariatePolynomial, - _b: &MultivariatePolynomial, - ) -> Option<( - MultivariatePolynomial, - MultivariatePolynomial, - MultivariatePolynomial, - )> { - None - } - - fn gcd( - a: &MultivariatePolynomial, - b: &MultivariatePolynomial, - vars: &[usize], - bounds: &mut [E], - tight_bounds: &mut [E], - ) -> MultivariatePolynomial { - assert!(!a.is_zero() || !b.is_zero()); - MultivariatePolynomial::gcd_shape_modular(a, b, vars, bounds, tight_bounds).unwrap() - } - - fn get_gcd_var_bounds( - a: &MultivariatePolynomial, - b: &MultivariatePolynomial, - vars: &[usize], - loose_bounds: &[E], - ) -> SmallVec<[E; INLINED_EXPONENTS]> { - let mut tight_bounds: SmallVec<[_; INLINED_EXPONENTS]> = loose_bounds.into(); - for var in vars { - let vvars: SmallVec<[usize; INLINED_EXPONENTS]> = - vars.iter().filter(|i| *i != var).cloned().collect(); - tight_bounds[*var] = MultivariatePolynomial::get_gcd_var_bound(a, b, &vvars, *var); - } - tight_bounds - } - - fn gcd_multiple(f: Vec>) -> MultivariatePolynomial { - MultivariatePolynomial::repeated_gcd(f) - } - - fn normalize(a: MultivariatePolynomial) -> MultivariatePolynomial { - a.make_monic() - } -} diff --git a/src/poly/polynomial.rs b/src/poly/polynomial.rs index 84f09e0..93290e7 100755 --- a/src/poly/polynomial.rs +++ b/src/poly/polynomial.rs @@ -1790,7 +1790,7 @@ impl MultivariatePolynomial { c2.exponents = c .exponents_iter() - .map(|x| x[var_index].to_u32() as u8) + .map(|x| x[var_index].to_u32() as u16) .collect(); c2.coefficients = c.coefficients;