From 24bb31e2e31492f9ff8b048c541435aecef287a6 Mon Sep 17 00:00:00 2001 From: Ben Ruijl Date: Wed, 17 Jul 2024 13:27:04 +0200 Subject: [PATCH] Add Fraction fields - Remove PartialOrd for ring elements and replace it with an internal ordering --- src/atom.rs | 59 +++ src/domains.rs | 49 +- src/domains/algebraic_number.rs | 10 +- src/domains/atom.rs | 8 +- src/domains/factorized_rational_polynomial.rs | 6 +- src/domains/finite_field.rs | 8 +- src/domains/integer.rs | 16 +- src/domains/rational.rs | 448 +++++++++++++++++- src/domains/rational_polynomial.rs | 38 +- src/poly/factor.rs | 35 +- src/poly/polynomial.rs | 35 +- src/poly/series.rs | 7 +- src/poly/univariate.rs | 11 +- 13 files changed, 647 insertions(+), 83 deletions(-) diff --git a/src/atom.rs b/src/atom.rs index 7c0e43d9..ccb2b4db 100644 --- a/src/atom.rs +++ b/src/atom.rs @@ -765,6 +765,20 @@ impl Atom { } } +impl std::ops::Add for Atom { + type Output = Atom; + + fn add(self, mut rhs: Atom) -> Atom { + Workspace::get_local().with(|ws| { + let mut t = ws.new_atom(); + self.as_view().add_with_ws_into(ws, rhs.as_view(), &mut t); + std::mem::swap(&mut rhs, &mut t); + }); + + rhs + } +} + impl std::ops::Add for &Atom { type Output = Atom; @@ -779,6 +793,23 @@ impl std::ops::Add for &Atom { } } +impl std::ops::Sub for Atom { + type Output = Atom; + + fn sub(self, mut rhs: Atom) -> Atom { + Workspace::get_local().with(|ws| { + let mut t = ws.new_atom(); + self.as_view() + .sub_no_norm(ws, rhs.as_view()) + .as_view() + .normalize(ws, &mut t); + std::mem::swap(&mut rhs, &mut t); + }); + + rhs + } +} + impl std::ops::Sub for &Atom { type Output = Atom; @@ -810,6 +841,34 @@ impl std::ops::Mul for &Atom { } } +impl std::ops::Mul for Atom { + type Output = Atom; + + fn mul(self, mut rhs: Atom) -> Atom { + Workspace::get_local().with(|ws| { + let mut t = ws.new_atom(); + self.as_view().mul_with_ws_into(ws, rhs.as_view(), &mut t); + std::mem::swap(&mut rhs, &mut t); + }); + + rhs + } +} + +impl std::ops::Div for Atom { + type Output = Atom; + + fn div(self, mut rhs: Atom) -> Atom { + Workspace::get_local().with(|ws| { + let mut t = ws.new_atom(); + self.as_view().div_with_ws_into(ws, rhs.as_view(), &mut t); + std::mem::swap(&mut rhs, &mut t); + }); + + rhs + } +} + impl std::ops::Div for &Atom { type Output = Atom; diff --git a/src/domains.rs b/src/domains.rs index f21c186d..76e42341 100644 --- a/src/domains.rs +++ b/src/domains.rs @@ -14,8 +14,55 @@ use integer::Integer; use crate::printer::PrintOptions; +pub trait InternalOrdering { + /// Compare two elements using an internal ordering. + fn internal_cmp(&self, other: &Self) -> std::cmp::Ordering; +} + +macro_rules! impl_internal_ordering { + ($($t:ty),*) => { + $( + impl InternalOrdering for $t { + fn internal_cmp(&self, other: &Self) -> std::cmp::Ordering { + self.cmp(other) + } + } + )* + }; +} + +impl_internal_ordering!(u8); +impl_internal_ordering!(u64); + +macro_rules! impl_internal_ordering_range { + ($($t:ty),*) => { + $( + impl InternalOrdering for $t { + fn internal_cmp(&self, other: &Self) -> std::cmp::Ordering { + match self.len().cmp(&other.len()) { + std::cmp::Ordering::Equal => (), + ord => return ord, + } + + for (i, j) in self.iter().zip(other) { + match i.internal_cmp(&j) { + std::cmp::Ordering::Equal => {} + ord => return ord, + } + } + + std::cmp::Ordering::Equal + } + } + )* + }; +} + +impl_internal_ordering_range!([T]); +impl_internal_ordering_range!(Vec); + pub trait Ring: Clone + PartialEq + Eq + Hash + Debug + Display { - type Element: Clone + PartialEq + Eq + Hash + PartialOrd + Debug; + type Element: Clone + PartialEq + Eq + Hash + InternalOrdering + Debug; fn add(&self, a: &Self::Element, b: &Self::Element) -> Self::Element; fn sub(&self, a: &Self::Element, b: &Self::Element) -> Self::Element; diff --git a/src/domains/algebraic_number.rs b/src/domains/algebraic_number.rs index 3009cc4e..3fc43e41 100644 --- a/src/domains/algebraic_number.rs +++ b/src/domains/algebraic_number.rs @@ -18,12 +18,12 @@ use super::{ }, integer::Integer, rational::Rational, - EuclideanDomain, Field, Ring, + EuclideanDomain, Field, InternalOrdering, Ring, }; /// An algebraic number ring, with a monic, irreducible defining polynomial. // TODO: make special case for degree two and three and hardcode the multiplication table -#[derive(Clone, PartialEq, Eq, PartialOrd, Hash)] +#[derive(Clone, PartialEq, Eq, Hash)] pub struct AlgebraicExtension { poly: Arc>, // TODO: convert to univariate polynomial } @@ -334,9 +334,9 @@ pub struct AlgebraicNumber { pub(crate) poly: MultivariatePolynomial, } -impl PartialOrd for AlgebraicNumber { - fn partial_cmp(&self, other: &Self) -> Option { - self.poly.partial_cmp(&other.poly) +impl InternalOrdering for AlgebraicNumber { + fn internal_cmp(&self, other: &Self) -> std::cmp::Ordering { + self.poly.internal_cmp(&other.poly) } } diff --git a/src/domains/atom.rs b/src/domains/atom.rs index bcf99039..1dddc8a4 100644 --- a/src/domains/atom.rs +++ b/src/domains/atom.rs @@ -1,6 +1,6 @@ use crate::atom::{Atom, AtomView}; -use super::{integer::Integer, EuclideanDomain, Field, Ring}; +use super::{integer::Integer, EuclideanDomain, Field, InternalOrdering, Ring}; use rand::Rng; @@ -31,6 +31,12 @@ impl std::fmt::Debug for AtomField { } } +impl InternalOrdering for Atom { + fn internal_cmp(&self, other: &Self) -> std::cmp::Ordering { + self.cmp(other) + } +} + impl Ring for AtomField { type Element = Atom; diff --git a/src/domains/factorized_rational_polynomial.rs b/src/domains/factorized_rational_polynomial.rs index 41701589..9ec7518a 100644 --- a/src/domains/factorized_rational_polynomial.rs +++ b/src/domains/factorized_rational_polynomial.rs @@ -19,7 +19,7 @@ use super::{ finite_field::{FiniteField, FiniteFieldCore, FiniteFieldWorkspace, ToFiniteField}, integer::{Integer, IntegerRing, Z}, rational::RationalField, - EuclideanDomain, Field, Ring, + EuclideanDomain, Field, InternalOrdering, Ring, }; #[derive(Clone, PartialEq, Eq, Hash, Debug)] @@ -71,9 +71,9 @@ pub struct FactorizedRationalPolynomial { pub denominators: Vec<(MultivariatePolynomial, usize)>, // TODO: sort factors? } -impl PartialOrd for FactorizedRationalPolynomial { +impl InternalOrdering for FactorizedRationalPolynomial { /// An ordering of rational polynomials that has no intuitive meaning. - fn partial_cmp(&self, _other: &Self) -> Option { + fn internal_cmp(&self, _other: &Self) -> Ordering { todo!() } } diff --git a/src/domains/finite_field.rs b/src/domains/finite_field.rs index 750192fd..f0bccb93 100644 --- a/src/domains/finite_field.rs +++ b/src/domains/finite_field.rs @@ -10,7 +10,7 @@ use crate::printer::PrintOptions; use super::algebraic_number::AlgebraicExtension; use super::integer::Z; -use super::{EuclideanDomain, Field, Ring}; +use super::{EuclideanDomain, Field, InternalOrdering, Ring}; const HENSEL_LIFTING_MASK: [u8; 128] = [ 255, 85, 51, 73, 199, 93, 59, 17, 15, 229, 195, 89, 215, 237, 203, 33, 31, 117, 83, 105, 231, @@ -145,6 +145,12 @@ where #[derive(Debug, Copy, Clone, Hash, PartialEq, PartialOrd, Eq)] pub struct FiniteFieldElement(pub(crate) UField); +impl InternalOrdering for FiniteFieldElement { + fn internal_cmp(&self, other: &Self) -> std::cmp::Ordering { + self.partial_cmp(other).unwrap_or(std::cmp::Ordering::Equal) + } +} + pub trait FiniteFieldWorkspace: Clone + Display + Eq + Hash { /// Convert to u64. fn to_u64(&self) -> u64; diff --git a/src/domains/integer.rs b/src/domains/integer.rs index 9743f2ab..f7dc1a73 100644 --- a/src/domains/integer.rs +++ b/src/domains/integer.rs @@ -20,7 +20,7 @@ use super::{ Zp64, Z2, }, rational::Rational, - EuclideanDomain, Ring, + EuclideanDomain, InternalOrdering, Ring, }; pub const SMALL_PRIMES: [i64; 100] = [ @@ -59,6 +59,12 @@ pub enum Integer { Large(MultiPrecisionInteger), } +impl InternalOrdering for Integer { + fn internal_cmp(&self, other: &Self) -> std::cmp::Ordering { + Ord::cmp(self, &other) + } +} + macro_rules! from_with_cast { ($base: ty) => { impl From<$base> for Integer { @@ -478,7 +484,7 @@ impl Integer { .as_abs() .partial_cmp(&n2.unsigned_abs()) .unwrap_or(Ordering::Equal), - (_, _) => self.abs().cmp(&other.abs()), + (_, _) => Ord::cmp(&self.abs(), &other.abs()), } } @@ -1871,6 +1877,12 @@ impl MultiPrecisionIntegerRing { } } +impl InternalOrdering for MultiPrecisionInteger { + fn internal_cmp(&self, other: &Self) -> std::cmp::Ordering { + Ord::cmp(self, &other) + } +} + impl Ring for MultiPrecisionIntegerRing { type Element = MultiPrecisionInteger; diff --git a/src/domains/rational.rs b/src/domains/rational.rs index 597f2324..9970e68e 100644 --- a/src/domains/rational.rs +++ b/src/domains/rational.rs @@ -1,4 +1,5 @@ use std::{ + borrow::Cow, fmt::{Display, Error, Formatter, Write}, ops::{Add, AddAssign, Div, DivAssign, Mul, MulAssign, Neg, Sub, SubAssign}, }; @@ -9,14 +10,18 @@ use rug::{ Rational as MultiPrecisionRational, }; -use crate::{poly::gcd::LARGE_U32_PRIMES, printer::PrintOptions, utils}; +use crate::{ + poly::{gcd::LARGE_U32_PRIMES, polynomial::PolynomialRing, Exponent}, + printer::PrintOptions, + utils, +}; use super::{ finite_field::{ FiniteField, FiniteFieldCore, FiniteFieldWorkspace, ToFiniteField, Two, Zp, Z2, }, integer::{Integer, Z}, - EuclideanDomain, Field, Ring, + EuclideanDomain, Field, InternalOrdering, Ring, }; /// The field of rational numbers. @@ -40,6 +45,391 @@ impl RationalField { } } +/// The fraction field of `R`. +#[derive(Clone, PartialEq, Eq, Hash, Debug)] +pub struct FractionField { + ring: R, +} + +impl FractionField { + pub const fn new(ring: R) -> FractionField { + FractionField { ring } + } +} + +impl FractionField { + pub fn to_element_numerator(&self, numerator: R::Element) -> ::Element { + Fraction { + numerator, + denominator: self.ring.one(), + } + } +} + +impl FractionField { + pub fn to_element( + &self, + mut numerator: R::Element, + mut denominator: R::Element, + do_gcd: bool, + ) -> ::Element { + if do_gcd { + let g = self.ring.gcd(&numerator, &denominator); + if !self.ring.is_one(&g) { + numerator = self.ring.quot_rem(&numerator, &g).0; + denominator = self.ring.quot_rem(&denominator, &g).0; + } + } + + let f = self.ring.get_normalization_factor(&denominator); + + if self.ring.is_one(&f) { + Fraction { + numerator, + denominator, + } + } else { + Fraction { + numerator: self.ring.mul(&numerator, &f), + denominator: self.ring.mul(&denominator, &f), + } + } + } +} + +impl Display for FractionField { + fn fmt(&self, _f: &mut Formatter<'_>) -> std::fmt::Result { + Ok(()) + } +} + +pub trait FractionNormalization: Ring { + /// Get the factor that normalizes the element `a`. + /// - For a field, this is the inverse of `a`. + /// - For the integers, this is the sign of `a`. + /// - For a polynomial ring, this is the normalization factor of the leading coefficient. + fn get_normalization_factor(&self, a: &Self::Element) -> Self::Element; +} + +impl FractionNormalization for Z { + fn get_normalization_factor(&self, a: &Integer) -> Integer { + if *a < 0 { + (-1).into() + } else { + 1.into() + } + } +} + +impl FractionNormalization for PolynomialRing { + fn get_normalization_factor(&self, a: &Self::Element) -> Self::Element { + a.constant(a.field.get_normalization_factor(&a.lcoeff())) + } +} + +impl FractionNormalization for T { + fn get_normalization_factor(&self, a: &Self::Element) -> Self::Element { + self.inv(a) + } +} + +#[derive(Clone, PartialEq, Eq, Hash, Debug)] +pub struct Fraction { + numerator: R::Element, + denominator: R::Element, +} + +impl Fraction { + pub fn numerator(&self) -> &R::Element { + &self.numerator + } + + pub fn denominator(&self) -> &R::Element { + &self.denominator + } +} + +impl InternalOrdering for Fraction { + fn internal_cmp(&self, other: &Self) -> std::cmp::Ordering { + self.numerator + .internal_cmp(&other.numerator) + .then_with(|| self.denominator.internal_cmp(&other.denominator)) + } +} + +impl Ring for FractionField { + type Element = Fraction; + + fn add(&self, a: &Self::Element, b: &Self::Element) -> Self::Element { + let r = &self.ring; + + if a.denominator == b.denominator { + let num = r.add(&a.numerator, &b.numerator); + let g = r.gcd(&num, &a.denominator); + if !r.is_one(&g) { + return Fraction { + numerator: r.quot_rem(&num, &g).0, + denominator: a.denominator.clone(), + }; + } else { + return Fraction { + numerator: num, + denominator: a.denominator.clone(), + }; + } + } + + let denom_gcd = r.gcd(&a.denominator, &b.denominator); + + let mut a_den_red = Cow::Borrowed(&a.denominator); + let mut b_den_red = Cow::Borrowed(&b.denominator); + + if !r.is_one(&denom_gcd) { + a_den_red = Cow::Owned(r.quot_rem(&a.denominator, &denom_gcd).0); + b_den_red = Cow::Owned(r.quot_rem(&b.denominator, &denom_gcd).0); + } + + let num1 = r.mul(&a.numerator, &b_den_red); + let num2 = r.mul(&b.numerator, &a_den_red); + let mut num = r.add(&num1, &num2); + + // TODO: prefer small * large over medium * medium sized operations + // a_denom_red.as_ref() * &other.denominator may be faster + // TODO: add size hint trait with default implementation? + let mut den = r.mul(b_den_red.as_ref(), &a.denominator); + + let g = r.gcd(&num, &denom_gcd); + + if !r.is_one(&g) { + num = r.quot_rem(&num, &g).0; + den = r.quot_rem(&den, &g).0; + } + + Fraction { + numerator: num, + denominator: den, + } + } + + fn sub(&self, a: &Self::Element, b: &Self::Element) -> Self::Element { + // TODO: optimize + self.add(a, &self.neg(b)) + } + + fn mul(&self, a: &Self::Element, b: &Self::Element) -> Self::Element { + let r = &self.ring; + let gcd1 = r.gcd(&a.numerator, &b.denominator); + let gcd2 = r.gcd(&a.denominator, &b.numerator); + + if r.is_one(&gcd1) { + if r.is_one(&gcd2) { + Fraction { + numerator: r.mul(&a.numerator, &b.numerator), + denominator: r.mul(&a.denominator, &b.denominator), + } + } else { + Fraction { + numerator: r.mul(&a.numerator, &r.quot_rem(&b.numerator, &gcd2).0), + denominator: self + .ring + .mul(&r.quot_rem(&a.denominator, &gcd2).0, &b.denominator), + } + } + } else if r.is_one(&gcd2) { + Fraction { + numerator: r.mul(&r.quot_rem(&a.numerator, &gcd1).0, &b.numerator), + denominator: r.mul(&a.denominator, &r.quot_rem(&b.denominator, &gcd1).0), + } + } else { + Fraction { + numerator: r.mul( + &r.quot_rem(&a.numerator, &gcd1).0, + &r.quot_rem(&b.numerator, &gcd2).0, + ), + denominator: r.mul( + &r.quot_rem(&a.denominator, &gcd2).0, + &r.quot_rem(&b.denominator, &gcd1).0, + ), + } + } + } + + fn add_assign(&self, a: &mut Self::Element, b: &Self::Element) { + // TODO: optimize + *a = self.add(a, b); + } + + fn sub_assign(&self, a: &mut Self::Element, b: &Self::Element) { + *a = self.sub(a, b); + } + + 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)); + } + + fn neg(&self, a: &Self::Element) -> Self::Element { + Fraction { + numerator: self.ring.neg(&a.numerator), + denominator: a.denominator.clone(), + } + } + + fn zero(&self) -> Self::Element { + Fraction { + numerator: self.ring.zero(), + denominator: self.ring.one(), + } + } + + fn one(&self) -> Self::Element { + Fraction { + numerator: self.ring.one(), + denominator: self.ring.one(), + } + } + + #[inline] + fn nth(&self, n: u64) -> Self::Element { + Fraction { + numerator: self.ring.nth(n), + denominator: self.ring.one(), + } + } + + fn pow(&self, b: &Self::Element, e: u64) -> Self::Element { + Fraction { + numerator: self.ring.pow(&b.numerator, e), + denominator: self.ring.pow(&b.denominator, e), + } + } + + fn is_zero(a: &Self::Element) -> bool { + R::is_zero(&a.numerator) + } + + fn is_one(&self, a: &Self::Element) -> bool { + self.ring.is_one(&a.numerator) && self.ring.is_one(&a.denominator) + } + + fn one_is_gcd_unit() -> bool { + false + } + + fn characteristic(&self) -> Integer { + self.ring.characteristic() + } + + fn size(&self) -> Integer { + // TODO: this is an overestimate + self.ring.size() * self.ring.size() + } + + fn sample(&self, rng: &mut impl rand::RngCore, range: (i64, i64)) -> Self::Element { + Fraction { + numerator: self.ring.sample(rng, range), + denominator: self.ring.one(), + } + } + + fn fmt_display( + &self, + element: &Self::Element, + opts: &PrintOptions, + _in_product: bool, + f: &mut Formatter<'_>, + ) -> Result<(), Error> { + self.ring.fmt_display(&element.numerator, opts, true, f)?; + if !self.ring.is_one(&element.denominator) { + f.write_char('/')?; + self.ring.fmt_display(&element.denominator, opts, true, f)?; + } + + Ok(()) + } +} + +impl EuclideanDomain for FractionField { + fn rem(&self, _: &Self::Element, _: &Self::Element) -> Self::Element { + self.zero() + } + + fn quot_rem(&self, a: &Self::Element, b: &Self::Element) -> (Self::Element, Self::Element) { + (self.div(a, b), self.zero()) + } + + fn gcd(&self, a: &Self::Element, b: &Self::Element) -> Self::Element { + let gcd_num = self.ring.gcd(&a.numerator, &b.denominator); + let gcd_den = self.ring.gcd(&a.denominator, &b.numerator); + + let d1 = self.ring.quot_rem(&a.denominator, &gcd_den).0; + let lcm = self.ring.mul(&d1, &b.denominator); + + Fraction { + numerator: gcd_num, + denominator: lcm, + } + } +} + +impl Field for FractionField { + fn div(&self, a: &Self::Element, b: &Self::Element) -> Self::Element { + // TODO: optimize + self.mul(a, &self.inv(b)) + } + + fn div_assign(&self, a: &mut Self::Element, b: &Self::Element) { + *a = self.div(a, b); + } + + fn inv(&self, a: &Self::Element) -> Self::Element { + let f = self.ring.get_normalization_factor(&a.numerator); + + Fraction { + numerator: self.ring.mul(&a.denominator, &f), + denominator: self.ring.mul(&a.numerator, &f), + } + } +} + +impl PolynomialRing, E> { + pub fn to_rational_polynomial( + &self, + e: &::Element, + ) -> Fraction> { + let mut lcm = self.ring.ring.one(); + for x in &e.coefficients { + let g = self.ring.ring.gcd(&lcm, x.denominator()); + lcm = self + .ring + .ring + .mul(&lcm, &self.ring.ring.quot_rem(x.denominator(), &g).0); + } + + let e2 = e.map_coeff( + |c| { + self.ring.ring.mul( + &c.numerator, + &self.ring.ring.quot_rem(&lcm, &c.denominator).0, + ) + }, + self.ring.ring.clone(), + ); + + Fraction { + denominator: e2.constant(lcm), + numerator: e2, + } + } +} + /// A rational number. /// /// Explicit construction of `Rational::Natural` @@ -54,6 +444,12 @@ pub enum Rational { Large(MultiPrecisionRational), } +impl InternalOrdering for Rational { + fn internal_cmp(&self, other: &Self) -> std::cmp::Ordering { + Ord::cmp(self, &other) + } +} + impl From for Rational { #[inline] fn from(value: i32) -> Self { @@ -1056,7 +1452,15 @@ impl<'a> std::iter::Sum<&'a Self> for Rational { #[cfg(test)] mod test { - use crate::domains::rational::Rational; + use crate::{ + atom::Atom, + domains::{ + integer::Z, + rational::{FractionField, Rational, Q}, + Field, Ring, + }, + poly::polynomial::PolynomialRing, + }; #[test] fn rounding() { @@ -1084,4 +1488,42 @@ mod test { let res = r.round(&Rational::new(1, 100000000)); assert_eq!(res, (93343, 29712).into()); } + + #[test] + fn fraction_int() { + let f = FractionField::new(Z); + let b = f.neg(&f.nth(3)); + let d = f.div(&f.add(&f.nth(100), &b), &b); + assert_eq!(d, f.to_element((-97).into(), 3.into(), false)); + } + + #[test] + fn fraction_poly() { + let poly = Atom::parse("-3/2*x^2+1/5x+4") + .unwrap() + .to_polynomial::<_, u8>(&Q, None); + + let f = FractionField::new(Z); + let poly2 = poly.map_coeff( + |c| f.to_element(c.numerator(), c.denominator(), false), + f.clone(), + ); + + let p = PolynomialRing::new_from_poly(&poly2); + let rat = p.to_rational_polynomial(&poly2); + let f = FractionField::new(PolynomialRing::new_from_poly(&rat.numerator)); + + let b = f.neg(&f.nth(3)); + let c = f.add(&rat, &b); + let d = f.div(&c, &rat); + + let num = Atom::parse("-10-2*x+15*x^2") + .unwrap() + .to_polynomial::<_, u8>(&Z, None); + let den = Atom::parse("-40-2*x+15*x^2") + .unwrap() + .to_polynomial::<_, u8>(&Z, None); + + assert_eq!(d, f.to_element(num, den, false)); + } } diff --git a/src/domains/rational_polynomial.rs b/src/domains/rational_polynomial.rs index 121f6981..b4c1d15c 100644 --- a/src/domains/rational_polynomial.rs +++ b/src/domains/rational_polynomial.rs @@ -21,7 +21,7 @@ use super::{ finite_field::{FiniteField, FiniteFieldCore, FiniteFieldWorkspace, ToFiniteField}, integer::{Integer, IntegerRing, Z}, rational::RationalField, - EuclideanDomain, Field, Ring, + EuclideanDomain, Field, InternalOrdering, Ring, }; #[derive(Clone, PartialEq, Eq, Hash, Debug)] @@ -64,27 +64,23 @@ pub struct RationalPolynomial { pub denominator: MultivariatePolynomial, } -impl PartialOrd for RationalPolynomial { +impl InternalOrdering for RationalPolynomial { /// An ordering of rational polynomials that has no intuitive meaning. - fn partial_cmp(&self, other: &Self) -> Option { - Some( - self.numerator - .exponents - .cmp(&other.numerator.exponents) - .then_with(|| self.denominator.exponents.cmp(&other.denominator.exponents)) - .then_with(|| { - self.numerator - .coefficients - .partial_cmp(&other.numerator.coefficients) - .unwrap_or(Ordering::Equal) - }) - .then_with(|| { - self.denominator - .coefficients - .partial_cmp(&other.denominator.coefficients) - .unwrap_or(Ordering::Equal) - }), - ) + fn internal_cmp(&self, other: &Self) -> Ordering { + self.numerator + .exponents + .cmp(&other.numerator.exponents) + .then_with(|| self.denominator.exponents.cmp(&other.denominator.exponents)) + .then_with(|| { + self.numerator + .coefficients + .internal_cmp(&other.numerator.coefficients) + }) + .then_with(|| { + self.denominator + .coefficients + .internal_cmp(&other.denominator.coefficients) + }) } } diff --git a/src/poly/factor.rs b/src/poly/factor.rs index 088b0f94..feeca181 100644 --- a/src/poly/factor.rs +++ b/src/poly/factor.rs @@ -14,7 +14,7 @@ use crate::{ }, integer::{Integer, IntegerRing, Z}, rational::{RationalField, Q}, - EuclideanDomain, Field, Ring, + EuclideanDomain, Field, InternalOrdering, Ring, }, poly::gcd::LARGE_U32_PRIMES, utils, @@ -1165,7 +1165,7 @@ where univariate_factors.sort_by(|(_, _, a), (_, _, b)| { a.exponents .cmp(&b.exponents) - .then(a.coefficients.partial_cmp(&b.coefficients).unwrap()) + .then(a.coefficients.internal_cmp(&b.coefficients)) }); univariate_factors @@ -3377,6 +3377,7 @@ mod test { finite_field::{Zp, Z2}, integer::Z, rational::Q, + InternalOrdering, }, poly::factor::Factorize, }; @@ -3402,9 +3403,9 @@ mod test { ) }) .collect::>(); - res.sort_by(|a, b| a.partial_cmp(&b).unwrap()); + res.sort_by(|a, b| a.0.internal_cmp(&b.0).then(a.1.cmp(&b.1))); let mut r = poly.square_free_factorization(); - r.sort_by(|a, b| a.partial_cmp(&b).unwrap()); + r.sort_by(|a, b| a.0.internal_cmp(&b.0).then(a.1.cmp(&b.1))); assert_eq!(r, res); } @@ -3431,9 +3432,9 @@ mod test { }) .collect::>(); - res.sort_by(|a, b| a.partial_cmp(&b).unwrap()); + res.sort_by(|a, b| a.0.internal_cmp(&b.0).then(a.1.cmp(&b.1))); let mut r = poly.factor(); - r.sort_by(|a, b| a.partial_cmp(&b).unwrap()); + r.sort_by(|a, b| a.0.internal_cmp(&b.0).then(a.1.cmp(&b.1))); assert_eq!(r, res); } @@ -3464,9 +3465,9 @@ mod test { }) .collect::>(); - res.sort_by(|a, b| a.partial_cmp(&b).unwrap()); + res.sort_by(|a, b| a.0.internal_cmp(&b.0).then(a.1.cmp(&b.1))); let mut r = poly.square_free_factorization(); - r.sort_by(|a, b| a.partial_cmp(&b).unwrap()); + r.sort_by(|a, b| a.0.internal_cmp(&b.0).then(a.1.cmp(&b.1))); assert_eq!(r, res); } @@ -3497,9 +3498,9 @@ mod test { }) .collect::>(); - res.sort_by(|a, b| a.partial_cmp(&b).unwrap()); + res.sort_by(|a, b| a.0.internal_cmp(&b.0).then(a.1.cmp(&b.1))); let mut r = poly.factor(); - r.sort_by(|a, b| a.partial_cmp(&b).unwrap()); + r.sort_by(|a, b| a.0.internal_cmp(&b.0).then(a.1.cmp(&b.1))); assert_eq!(r, res); } @@ -3537,9 +3538,9 @@ mod test { }) .collect::>(); - res.sort_by(|a, b| a.partial_cmp(&b).unwrap()); + res.sort_by(|a, b| a.0.internal_cmp(&b.0).then(a.1.cmp(&b.1))); let mut r = poly.factor(); - r.sort_by(|a, b| a.partial_cmp(&b).unwrap()); + r.sort_by(|a, b| a.0.internal_cmp(&b.0).then(a.1.cmp(&b.1))); assert_eq!(r, res); } @@ -3568,9 +3569,9 @@ mod test { }) .collect::>(); - res.sort_by(|a, b| a.partial_cmp(&b).unwrap()); + res.sort_by(|a, b| a.0.internal_cmp(&b.0).then(a.1.cmp(&b.1))); let mut r = poly.factor(); - r.sort_by(|a, b| a.partial_cmp(&b).unwrap()); + r.sort_by(|a, b| a.0.internal_cmp(&b.0).then(a.1.cmp(&b.1))); assert_eq!(r, res); } @@ -3598,9 +3599,9 @@ mod test { }) .collect::>(); - res.sort_by(|a, b| a.partial_cmp(&b).unwrap()); + res.sort_by(|a, b| a.0.internal_cmp(&b.0).then(a.1.cmp(&b.1))); let mut r = poly.factor(); - r.sort_by(|a, b| a.partial_cmp(&b).unwrap()); + r.sort_by(|a, b| a.0.internal_cmp(&b.0).then(a.1.cmp(&b.1))); assert_eq!(r, res); } @@ -3636,7 +3637,7 @@ mod test { .to_polynomial::<_, u8>(&Q, a.get_vars().clone().into()) .to_number_field(&f); - factors.sort_by(|a, b| a.partial_cmp(&b).unwrap()); + factors.sort_by(|a, b| a.0.internal_cmp(&b.0).then(a.1.cmp(&b.1))); assert_eq!(factors, vec![(f1, 1), (f2, 1)]) } diff --git a/src/poly/polynomial.rs b/src/poly/polynomial.rs index 4b864ad4..dd1bf6cb 100755 --- a/src/poly/polynomial.rs +++ b/src/poly/polynomial.rs @@ -11,7 +11,7 @@ use std::sync::Arc; use crate::domains::algebraic_number::AlgebraicExtension; use crate::domains::integer::{Integer, IntegerRing}; use crate::domains::rational::RationalField; -use crate::domains::{EuclideanDomain, Field, Ring}; +use crate::domains::{EuclideanDomain, Field, InternalOrdering, Ring}; use crate::printer::{PolynomialPrinter, PrintOptions}; use super::gcd::PolynomialGCD; @@ -484,11 +484,7 @@ impl MultivariatePolynomial { self.exponents = newexp; // reconstruct 'other' with correct monomial ordering - let mut newother = Self::new( - &other.field, - other.nterms().into(), - self.variables.clone(), - ); + let mut newother = Self::new(&other.field, other.nterms().into(), self.variables.clone()); let mut newexp = vec![E::zero(); self.nvars()]; for t in other.into_iter() { for c in &mut newexp { @@ -776,14 +772,11 @@ impl std::hash::Hash for MultivariatePol impl Eq for MultivariatePolynomial {} -impl PartialOrd for MultivariatePolynomial { +impl InternalOrdering for MultivariatePolynomial { /// An ordering of polynomials that has no intuitive meaning. - fn partial_cmp(&self, other: &Self) -> Option { - Some(self.exponents.cmp(&other.exponents).then_with(|| { - self.coefficients - .partial_cmp(&other.coefficients) - .unwrap_or(Ordering::Equal) - })) + fn internal_cmp(&self, other: &Self) -> Ordering { + Ord::cmp(&self.exponents, &other.exponents) + .then_with(|| self.coefficients.internal_cmp(&other.coefficients)) } } @@ -1432,11 +1425,8 @@ impl MultivariatePolynomial { let mut indices: Vec = (0..self.nterms()).collect(); indices.sort_unstable_by_key(|&i| &new_exp[i * order.len()..(i + 1) * order.len()]); - let mut res = MultivariatePolynomial::new( - &self.field, - self.nterms().into(), - self.variables.clone(), - ); + let mut res = + MultivariatePolynomial::new(&self.field, self.nterms().into(), self.variables.clone()); for i in indices { res.append_monomial( @@ -2828,7 +2818,7 @@ impl MultivariatePolynomial { /// Each exponent is limited to 32767 if there are 5 or fewer variables, /// or 127 if there are 8 or fewer variables, such that the last bit per byte can /// be used to check for subtraction overflow, serving as a division test. - pub fn heap_division_packed_exp( + fn heap_division_packed_exp( &self, div: &MultivariatePolynomial, abort_on_remainder: bool, @@ -3376,11 +3366,8 @@ impl MultivariatePolynomial, E> { let var_index = var_map.iter().position(|x| x == var).unwrap(); - let mut poly = MultivariatePolynomial::new( - &self.field.poly().field, - self.nterms().into(), - var_map, - ); + let mut poly = + MultivariatePolynomial::new(&self.field.poly().field, self.nterms().into(), var_map); let mut exp = vec![E::zero(); poly.nvars()]; for t in self { exp[..self.nvars()].copy_from_slice(t.exponents); diff --git a/src/poly/series.rs b/src/poly/series.rs index cd5d6cc6..5394a100 100644 --- a/src/poly/series.rs +++ b/src/poly/series.rs @@ -9,7 +9,8 @@ use crate::{ atom::{Atom, AtomView, FunctionBuilder}, coefficient::CoefficientView, domains::{ - atom::AtomField, integer::Integer, rational::Rational, EuclideanDomain, Ring, RingPrinter, + atom::AtomField, integer::Integer, rational::Rational, EuclideanDomain, InternalOrdering, + Ring, RingPrinter, }, printer::PrintOptions, state::State, @@ -553,9 +554,9 @@ impl std::hash::Hash for Series { impl Eq for Series {} impl PartialOrd for Series { - /// An ordering of seriess that has no intuitive meaning. + /// An ordering of series that has no intuitive meaning. fn partial_cmp(&self, other: &Self) -> Option { - self.coefficients.partial_cmp(&other.coefficients) + Some(self.coefficients.internal_cmp(&other.coefficients)) } } diff --git a/src/poly/univariate.rs b/src/poly/univariate.rs index 60f43667..4b40c0e4 100644 --- a/src/poly/univariate.rs +++ b/src/poly/univariate.rs @@ -5,7 +5,7 @@ use std::{ }; use crate::{ - domains::{integer::Integer, EuclideanDomain, Field, Ring, RingPrinter}, + domains::{integer::Integer, EuclideanDomain, Field, InternalOrdering, Ring, RingPrinter}, printer::PrintOptions, }; @@ -168,6 +168,13 @@ pub struct UnivariatePolynomial { pub field: F, } +impl InternalOrdering for UnivariatePolynomial { + /// An ordering of polynomials that has no intuitive meaning. + fn internal_cmp(&self, other: &Self) -> std::cmp::Ordering { + self.coefficients.internal_cmp(&other.coefficients) + } +} + impl std::fmt::Debug for UnivariatePolynomial { fn fmt(&self, f: &mut std::fmt::Formatter) -> std::fmt::Result { if self.is_zero() { @@ -544,7 +551,7 @@ impl Eq for UnivariatePolynomial {} impl PartialOrd for UnivariatePolynomial { /// An ordering of polynomials that has no intuitive meaning. fn partial_cmp(&self, other: &Self) -> Option { - self.coefficients.partial_cmp(&other.coefficients) + Some(self.coefficients.internal_cmp(&other.coefficients)) } }