diff --git a/examples/solve_linear_system.rs b/examples/solve_linear_system.rs index 7b752c5b..4ba3992a 100644 --- a/examples/solve_linear_system.rs +++ b/examples/solve_linear_system.rs @@ -58,7 +58,7 @@ fn solve_from_matrix() { }) .collect(); - let field = RationalPolynomialField::new_from_poly(&rhs_rat[0].numerator); + let field = RationalPolynomialField::from_poly(&rhs_rat[0].numerator); let m = Matrix::from_linear( system_rat, system.len() as u32, diff --git a/src/api/python.rs b/src/api/python.rs index 590bcc4a..999cd7ae 100644 --- a/src/api/python.rs +++ b/src/api/python.rs @@ -4368,17 +4368,8 @@ impl PythonPolynomial { /// Add two polynomials `self and `rhs`, returning the result. pub fn __add__(&self, rhs: Self) -> Self { - if self.poly.get_vars_ref() == rhs.poly.get_vars_ref() { - Self { - poly: self.poly.clone() + rhs.poly.clone(), - } - } else { - let mut new_self = self.poly.clone(); - let mut new_rhs = rhs.poly.clone(); - new_self.unify_variables(&mut new_rhs); - Self { - poly: new_self + new_rhs, - } + Self { + poly: self.poly.clone() + rhs.poly.clone(), } } @@ -4389,31 +4380,14 @@ impl PythonPolynomial { /// Multiply two polynomials `self and `rhs`, returning the result. pub fn __mul__(&self, rhs: Self) -> Self { - if self.poly.get_vars_ref() == rhs.poly.get_vars_ref() { - Self { - poly: &self.poly * &rhs.poly, - } - } else { - let mut new_self = self.poly.clone(); - let mut new_rhs = rhs.poly.clone(); - new_self.unify_variables(&mut new_rhs); - Self { - poly: new_self * &new_rhs, - } + Self { + poly: &self.poly * &rhs.poly, } } /// Divide the polynomial `self` by `rhs` if possible, returning the result. pub fn __truediv__(&self, rhs: Self) -> PyResult { - let (q, r) = if self.poly.get_vars_ref() == rhs.poly.get_vars_ref() { - self.poly.quot_rem(&rhs.poly, false) - } else { - let mut new_self = self.poly.clone(); - let mut new_rhs = rhs.poly.clone(); - new_self.unify_variables(&mut new_rhs); - - new_self.quot_rem(&new_rhs, false) - }; + let (q, r) = self.poly.quot_rem(&rhs.poly, false); if r.is_zero() { Ok(Self { poly: q }) @@ -4429,16 +4403,8 @@ impl PythonPolynomial { pub fn quot_rem(&self, rhs: Self) -> PyResult<(Self, Self)> { if rhs.poly.is_zero() { Err(exceptions::PyValueError::new_err("Division by zero")) - } else if self.poly.get_vars_ref() == rhs.poly.get_vars_ref() { - let (q, r) = self.poly.quot_rem(&rhs.poly, false); - Ok((Self { poly: q }, Self { poly: r })) } else { - let mut new_self = self.poly.clone(); - let mut new_rhs = rhs.poly.clone(); - new_self.unify_variables(&mut new_rhs); - - let (q, r) = new_self.quot_rem(&new_rhs, false); - + let (q, r) = self.poly.quot_rem(&rhs.poly, false); Ok((Self { poly: q }, Self { poly: r })) } } @@ -4454,34 +4420,17 @@ impl PythonPolynomial { pub fn __mod__(&self, rhs: Self) -> PyResult { if rhs.poly.is_zero() { Err(exceptions::PyValueError::new_err("Division by zero")) - } else if self.poly.get_vars_ref() == rhs.poly.get_vars_ref() { - Ok(Self { - poly: self.poly.rem(&rhs.poly), - }) } else { - let mut new_self = self.poly.clone(); - let mut new_rhs = rhs.poly.clone(); - new_self.unify_variables(&mut new_rhs); - Ok(Self { - poly: new_self.rem(&new_rhs), + poly: self.poly.rem(&rhs.poly), }) } } /// Compute the greatest common divisor (GCD) of two polynomials. pub fn gcd(&self, rhs: Self) -> Self { - if self.poly.get_vars_ref() == rhs.poly.get_vars_ref() { - Self { - poly: self.poly.gcd(&rhs.poly), - } - } else { - let mut new_self = self.poly.clone(); - let mut new_rhs = rhs.poly.clone(); - new_self.unify_variables(&mut new_rhs); - Self { - poly: new_self.gcd(&new_rhs), - } + Self { + poly: self.poly.gcd(&rhs.poly), } } @@ -5138,17 +5087,8 @@ impl PythonFiniteFieldPolynomial { /// Add two polynomials `self and `rhs`, returning the result. pub fn __add__(&self, rhs: Self) -> Self { - if self.poly.get_vars_ref() == rhs.poly.get_vars_ref() { - Self { - poly: self.poly.clone() + rhs.poly.clone(), - } - } else { - let mut new_self = self.poly.clone(); - let mut new_rhs = rhs.poly.clone(); - new_self.unify_variables(&mut new_rhs); - Self { - poly: new_self + new_rhs, - } + Self { + poly: self.poly.clone() + rhs.poly.clone(), } } @@ -5159,31 +5099,14 @@ impl PythonFiniteFieldPolynomial { /// Multiply two polynomials `self and `rhs`, returning the result. pub fn __mul__(&self, rhs: Self) -> Self { - if self.poly.get_vars_ref() == rhs.poly.get_vars_ref() { - Self { - poly: &self.poly * &rhs.poly, - } - } else { - let mut new_self = self.poly.clone(); - let mut new_rhs = rhs.poly.clone(); - new_self.unify_variables(&mut new_rhs); - Self { - poly: new_self * &new_rhs, - } + Self { + poly: &self.poly * &rhs.poly, } } /// Divide the polynomial `self` by `rhs` if possible, returning the result. pub fn __truediv__(&self, rhs: Self) -> PyResult { - let (q, r) = if self.poly.get_vars_ref() == rhs.poly.get_vars_ref() { - self.poly.quot_rem(&rhs.poly, false) - } else { - let mut new_self = self.poly.clone(); - let mut new_rhs = rhs.poly.clone(); - new_self.unify_variables(&mut new_rhs); - - new_self.quot_rem(&new_rhs, false) - }; + let (q, r) = self.poly.quot_rem(&rhs.poly, false); if r.is_zero() { Ok(Self { poly: q }) @@ -5199,16 +5122,8 @@ impl PythonFiniteFieldPolynomial { pub fn quot_rem(&self, rhs: Self) -> PyResult<(Self, Self)> { if rhs.poly.is_zero() { Err(exceptions::PyValueError::new_err("Division by zero")) - } else if self.poly.get_vars_ref() == rhs.poly.get_vars_ref() { - let (q, r) = self.poly.quot_rem(&rhs.poly, false); - Ok((Self { poly: q }, Self { poly: r })) } else { - let mut new_self = self.poly.clone(); - let mut new_rhs = rhs.poly.clone(); - new_self.unify_variables(&mut new_rhs); - - let (q, r) = new_self.quot_rem(&new_rhs, false); - + let (q, r) = self.poly.quot_rem(&rhs.poly, false); Ok((Self { poly: q }, Self { poly: r })) } } @@ -5224,34 +5139,17 @@ impl PythonFiniteFieldPolynomial { pub fn __mod__(&self, rhs: Self) -> PyResult { if rhs.poly.is_zero() { Err(exceptions::PyValueError::new_err("Division by zero")) - } else if self.poly.get_vars_ref() == rhs.poly.get_vars_ref() { - Ok(Self { - poly: self.poly.rem(&rhs.poly), - }) } else { - let mut new_self = self.poly.clone(); - let mut new_rhs = rhs.poly.clone(); - new_self.unify_variables(&mut new_rhs); - Ok(Self { - poly: new_self.rem(&new_rhs), + poly: self.poly.rem(&rhs.poly), }) } } /// Compute the greatest common divisor (GCD) of two polynomials. pub fn gcd(&self, rhs: Self) -> Self { - if self.poly.get_vars_ref() == rhs.poly.get_vars_ref() { - Self { - poly: self.poly.gcd(&rhs.poly), - } - } else { - let mut new_self = self.poly.clone(); - let mut new_rhs = rhs.poly.clone(); - new_self.unify_variables(&mut new_rhs); - Self { - poly: new_self.gcd(&new_rhs), - } + Self { + poly: self.poly.gcd(&rhs.poly), } } @@ -5608,7 +5506,7 @@ impl PythonFiniteFieldPolynomial { /// Convert the polynomial to an expression. pub fn to_expression(&self) -> PyResult { let p = self.poly.map_coeff( - |c| Integer::from_finite_field(&self.poly.field, c.clone()), + |c| Integer::from_finite_field(&self.poly.ring, c.clone()), IntegerRing::new(), ); @@ -5759,17 +5657,8 @@ impl PythonPrimeTwoPolynomial { /// Add two polynomials `self and `rhs`, returning the result. pub fn __add__(&self, rhs: Self) -> Self { - if self.poly.get_vars_ref() == rhs.poly.get_vars_ref() { - Self { - poly: self.poly.clone() + rhs.poly.clone(), - } - } else { - let mut new_self = self.poly.clone(); - let mut new_rhs = rhs.poly.clone(); - new_self.unify_variables(&mut new_rhs); - Self { - poly: new_self + new_rhs, - } + Self { + poly: self.poly.clone() + rhs.poly.clone(), } } @@ -5780,31 +5669,14 @@ impl PythonPrimeTwoPolynomial { /// Multiply two polynomials `self and `rhs`, returning the result. pub fn __mul__(&self, rhs: Self) -> Self { - if self.poly.get_vars_ref() == rhs.poly.get_vars_ref() { - Self { - poly: &self.poly * &rhs.poly, - } - } else { - let mut new_self = self.poly.clone(); - let mut new_rhs = rhs.poly.clone(); - new_self.unify_variables(&mut new_rhs); - Self { - poly: new_self * &new_rhs, - } + Self { + poly: &self.poly * &rhs.poly, } } /// Divide the polynomial `self` by `rhs` if possible, returning the result. pub fn __truediv__(&self, rhs: Self) -> PyResult { - let (q, r) = if self.poly.get_vars_ref() == rhs.poly.get_vars_ref() { - self.poly.quot_rem(&rhs.poly, false) - } else { - let mut new_self = self.poly.clone(); - let mut new_rhs = rhs.poly.clone(); - new_self.unify_variables(&mut new_rhs); - - new_self.quot_rem(&new_rhs, false) - }; + let (q, r) = self.poly.quot_rem(&rhs.poly, false); if r.is_zero() { Ok(Self { poly: q }) @@ -5820,16 +5692,8 @@ impl PythonPrimeTwoPolynomial { pub fn quot_rem(&self, rhs: Self) -> PyResult<(Self, Self)> { if rhs.poly.is_zero() { Err(exceptions::PyValueError::new_err("Division by zero")) - } else if self.poly.get_vars_ref() == rhs.poly.get_vars_ref() { - let (q, r) = self.poly.quot_rem(&rhs.poly, false); - Ok((Self { poly: q }, Self { poly: r })) } else { - let mut new_self = self.poly.clone(); - let mut new_rhs = rhs.poly.clone(); - new_self.unify_variables(&mut new_rhs); - - let (q, r) = new_self.quot_rem(&new_rhs, false); - + let (q, r) = self.poly.quot_rem(&rhs.poly, false); Ok((Self { poly: q }, Self { poly: r })) } } @@ -5845,34 +5709,17 @@ impl PythonPrimeTwoPolynomial { pub fn __mod__(&self, rhs: Self) -> PyResult { if rhs.poly.is_zero() { Err(exceptions::PyValueError::new_err("Division by zero")) - } else if self.poly.get_vars_ref() == rhs.poly.get_vars_ref() { - Ok(Self { - poly: self.poly.rem(&rhs.poly), - }) } else { - let mut new_self = self.poly.clone(); - let mut new_rhs = rhs.poly.clone(); - new_self.unify_variables(&mut new_rhs); - Ok(Self { - poly: new_self.rem(&new_rhs), + poly: self.poly.rem(&rhs.poly), }) } } /// Compute the greatest common divisor (GCD) of two polynomials. pub fn gcd(&self, rhs: Self) -> Self { - if self.poly.get_vars_ref() == rhs.poly.get_vars_ref() { - Self { - poly: self.poly.gcd(&rhs.poly), - } - } else { - let mut new_self = self.poly.clone(); - let mut new_rhs = rhs.poly.clone(); - new_self.unify_variables(&mut new_rhs); - Self { - poly: new_self.gcd(&new_rhs), - } + Self { + poly: self.poly.gcd(&rhs.poly), } } @@ -6345,17 +6192,8 @@ impl PythonGaloisFieldPrimeTwoPolynomial { /// Add two polynomials `self and `rhs`, returning the result. pub fn __add__(&self, rhs: Self) -> Self { - if self.poly.get_vars_ref() == rhs.poly.get_vars_ref() { - Self { - poly: self.poly.clone() + rhs.poly.clone(), - } - } else { - let mut new_self = self.poly.clone(); - let mut new_rhs = rhs.poly.clone(); - new_self.unify_variables(&mut new_rhs); - Self { - poly: new_self + new_rhs, - } + Self { + poly: self.poly.clone() + rhs.poly.clone(), } } @@ -6366,31 +6204,14 @@ impl PythonGaloisFieldPrimeTwoPolynomial { /// Multiply two polynomials `self and `rhs`, returning the result. pub fn __mul__(&self, rhs: Self) -> Self { - if self.poly.get_vars_ref() == rhs.poly.get_vars_ref() { - Self { - poly: &self.poly * &rhs.poly, - } - } else { - let mut new_self = self.poly.clone(); - let mut new_rhs = rhs.poly.clone(); - new_self.unify_variables(&mut new_rhs); - Self { - poly: new_self * &new_rhs, - } + Self { + poly: &self.poly * &rhs.poly, } } /// Divide the polynomial `self` by `rhs` if possible, returning the result. pub fn __truediv__(&self, rhs: Self) -> PyResult { - let (q, r) = if self.poly.get_vars_ref() == rhs.poly.get_vars_ref() { - self.poly.quot_rem(&rhs.poly, false) - } else { - let mut new_self = self.poly.clone(); - let mut new_rhs = rhs.poly.clone(); - new_self.unify_variables(&mut new_rhs); - - new_self.quot_rem(&new_rhs, false) - }; + let (q, r) = self.poly.quot_rem(&rhs.poly, false); if r.is_zero() { Ok(Self { poly: q }) @@ -6406,16 +6227,8 @@ impl PythonGaloisFieldPrimeTwoPolynomial { pub fn quot_rem(&self, rhs: Self) -> PyResult<(Self, Self)> { if rhs.poly.is_zero() { Err(exceptions::PyValueError::new_err("Division by zero")) - } else if self.poly.get_vars_ref() == rhs.poly.get_vars_ref() { - let (q, r) = self.poly.quot_rem(&rhs.poly, false); - Ok((Self { poly: q }, Self { poly: r })) } else { - let mut new_self = self.poly.clone(); - let mut new_rhs = rhs.poly.clone(); - new_self.unify_variables(&mut new_rhs); - - let (q, r) = new_self.quot_rem(&new_rhs, false); - + let (q, r) = self.poly.quot_rem(&rhs.poly, false); Ok((Self { poly: q }, Self { poly: r })) } } @@ -6431,34 +6244,17 @@ impl PythonGaloisFieldPrimeTwoPolynomial { pub fn __mod__(&self, rhs: Self) -> PyResult { if rhs.poly.is_zero() { Err(exceptions::PyValueError::new_err("Division by zero")) - } else if self.poly.get_vars_ref() == rhs.poly.get_vars_ref() { - Ok(Self { - poly: self.poly.rem(&rhs.poly), - }) } else { - let mut new_self = self.poly.clone(); - let mut new_rhs = rhs.poly.clone(); - new_self.unify_variables(&mut new_rhs); - Ok(Self { - poly: new_self.rem(&new_rhs), + poly: self.poly.rem(&rhs.poly), }) } } /// Compute the greatest common divisor (GCD) of two polynomials. pub fn gcd(&self, rhs: Self) -> Self { - if self.poly.get_vars_ref() == rhs.poly.get_vars_ref() { - Self { - poly: self.poly.gcd(&rhs.poly), - } - } else { - let mut new_self = self.poly.clone(); - let mut new_rhs = rhs.poly.clone(); - new_self.unify_variables(&mut new_rhs); - Self { - poly: new_self.gcd(&new_rhs), - } + Self { + poly: self.poly.gcd(&rhs.poly), } } @@ -6935,17 +6731,8 @@ impl PythonGaloisFieldPolynomial { /// Add two polynomials `self and `rhs`, returning the result. pub fn __add__(&self, rhs: Self) -> Self { - if self.poly.get_vars_ref() == rhs.poly.get_vars_ref() { - Self { - poly: self.poly.clone() + rhs.poly.clone(), - } - } else { - let mut new_self = self.poly.clone(); - let mut new_rhs = rhs.poly.clone(); - new_self.unify_variables(&mut new_rhs); - Self { - poly: new_self + new_rhs, - } + Self { + poly: self.poly.clone() + rhs.poly.clone(), } } @@ -6956,31 +6743,14 @@ impl PythonGaloisFieldPolynomial { /// Multiply two polynomials `self and `rhs`, returning the result. pub fn __mul__(&self, rhs: Self) -> Self { - if self.poly.get_vars_ref() == rhs.poly.get_vars_ref() { - Self { - poly: &self.poly * &rhs.poly, - } - } else { - let mut new_self = self.poly.clone(); - let mut new_rhs = rhs.poly.clone(); - new_self.unify_variables(&mut new_rhs); - Self { - poly: new_self * &new_rhs, - } + Self { + poly: &self.poly * &rhs.poly, } } /// Divide the polynomial `self` by `rhs` if possible, returning the result. pub fn __truediv__(&self, rhs: Self) -> PyResult { - let (q, r) = if self.poly.get_vars_ref() == rhs.poly.get_vars_ref() { - self.poly.quot_rem(&rhs.poly, false) - } else { - let mut new_self = self.poly.clone(); - let mut new_rhs = rhs.poly.clone(); - new_self.unify_variables(&mut new_rhs); - - new_self.quot_rem(&new_rhs, false) - }; + let (q, r) = self.poly.quot_rem(&rhs.poly, false); if r.is_zero() { Ok(Self { poly: q }) @@ -6996,16 +6766,8 @@ impl PythonGaloisFieldPolynomial { pub fn quot_rem(&self, rhs: Self) -> PyResult<(Self, Self)> { if rhs.poly.is_zero() { Err(exceptions::PyValueError::new_err("Division by zero")) - } else if self.poly.get_vars_ref() == rhs.poly.get_vars_ref() { - let (q, r) = self.poly.quot_rem(&rhs.poly, false); - Ok((Self { poly: q }, Self { poly: r })) } else { - let mut new_self = self.poly.clone(); - let mut new_rhs = rhs.poly.clone(); - new_self.unify_variables(&mut new_rhs); - - let (q, r) = new_self.quot_rem(&new_rhs, false); - + let (q, r) = self.poly.quot_rem(&rhs.poly, false); Ok((Self { poly: q }, Self { poly: r })) } } @@ -7021,34 +6783,17 @@ impl PythonGaloisFieldPolynomial { pub fn __mod__(&self, rhs: Self) -> PyResult { if rhs.poly.is_zero() { Err(exceptions::PyValueError::new_err("Division by zero")) - } else if self.poly.get_vars_ref() == rhs.poly.get_vars_ref() { - Ok(Self { - poly: self.poly.rem(&rhs.poly), - }) } else { - let mut new_self = self.poly.clone(); - let mut new_rhs = rhs.poly.clone(); - new_self.unify_variables(&mut new_rhs); - Ok(Self { - poly: new_self.rem(&new_rhs), + poly: self.poly.rem(&rhs.poly), }) } } /// Compute the greatest common divisor (GCD) of two polynomials. pub fn gcd(&self, rhs: Self) -> Self { - if self.poly.get_vars_ref() == rhs.poly.get_vars_ref() { - Self { - poly: self.poly.gcd(&rhs.poly), - } - } else { - let mut new_self = self.poly.clone(); - let mut new_rhs = rhs.poly.clone(); - new_self.unify_variables(&mut new_rhs); - Self { - poly: new_self.gcd(&new_rhs), - } + Self { + poly: self.poly.gcd(&rhs.poly), } } @@ -7374,7 +7119,7 @@ impl PythonGaloisFieldPolynomial { .poly .to_expression_with_coeff_map(|_, element, out| { let p = element.poly.map_coeff( - |c| Integer::from_finite_field(&element.poly.field, c.clone()), + |c| Integer::from_finite_field(&element.poly.ring, c.clone()), IntegerRing::new(), ); p.to_expression_into(out); @@ -7526,17 +7271,8 @@ impl PythonNumberFieldPolynomial { /// Add two polynomials `self and `rhs`, returning the result. pub fn __add__(&self, rhs: Self) -> Self { - if self.poly.get_vars_ref() == rhs.poly.get_vars_ref() { - Self { - poly: self.poly.clone() + rhs.poly.clone(), - } - } else { - let mut new_self = self.poly.clone(); - let mut new_rhs = rhs.poly.clone(); - new_self.unify_variables(&mut new_rhs); - Self { - poly: new_self + new_rhs, - } + Self { + poly: self.poly.clone() + rhs.poly.clone(), } } @@ -7547,31 +7283,14 @@ impl PythonNumberFieldPolynomial { /// Multiply two polynomials `self and `rhs`, returning the result. pub fn __mul__(&self, rhs: Self) -> Self { - if self.poly.get_vars_ref() == rhs.poly.get_vars_ref() { - Self { - poly: &self.poly * &rhs.poly, - } - } else { - let mut new_self = self.poly.clone(); - let mut new_rhs = rhs.poly.clone(); - new_self.unify_variables(&mut new_rhs); - Self { - poly: new_self * &new_rhs, - } + Self { + poly: &self.poly * &rhs.poly, } } /// Divide the polynomial `self` by `rhs` if possible, returning the result. pub fn __truediv__(&self, rhs: Self) -> PyResult { - let (q, r) = if self.poly.get_vars_ref() == rhs.poly.get_vars_ref() { - self.poly.quot_rem(&rhs.poly, false) - } else { - let mut new_self = self.poly.clone(); - let mut new_rhs = rhs.poly.clone(); - new_self.unify_variables(&mut new_rhs); - - new_self.quot_rem(&new_rhs, false) - }; + let (q, r) = self.poly.quot_rem(&rhs.poly, false); if r.is_zero() { Ok(Self { poly: q }) @@ -7587,16 +7306,8 @@ impl PythonNumberFieldPolynomial { pub fn quot_rem(&self, rhs: Self) -> PyResult<(Self, Self)> { if rhs.poly.is_zero() { Err(exceptions::PyValueError::new_err("Division by zero")) - } else if self.poly.get_vars_ref() == rhs.poly.get_vars_ref() { - let (q, r) = self.poly.quot_rem(&rhs.poly, false); - Ok((Self { poly: q }, Self { poly: r })) } else { - let mut new_self = self.poly.clone(); - let mut new_rhs = rhs.poly.clone(); - new_self.unify_variables(&mut new_rhs); - - let (q, r) = new_self.quot_rem(&new_rhs, false); - + let (q, r) = self.poly.quot_rem(&rhs.poly, false); Ok((Self { poly: q }, Self { poly: r })) } } @@ -7612,34 +7323,17 @@ impl PythonNumberFieldPolynomial { pub fn __mod__(&self, rhs: Self) -> PyResult { if rhs.poly.is_zero() { Err(exceptions::PyValueError::new_err("Division by zero")) - } else if self.poly.get_vars_ref() == rhs.poly.get_vars_ref() { - Ok(Self { - poly: self.poly.rem(&rhs.poly), - }) } else { - let mut new_self = self.poly.clone(); - let mut new_rhs = rhs.poly.clone(); - new_self.unify_variables(&mut new_rhs); - Ok(Self { - poly: new_self.rem(&new_rhs), + poly: self.poly.rem(&rhs.poly), }) } } /// Compute the greatest common divisor (GCD) of two polynomials. pub fn gcd(&self, rhs: Self) -> Self { - if self.poly.get_vars_ref() == rhs.poly.get_vars_ref() { - Self { - poly: self.poly.gcd(&rhs.poly), - } - } else { - let mut new_self = self.poly.clone(); - let mut new_rhs = rhs.poly.clone(); - new_self.unify_variables(&mut new_rhs); - Self { - poly: new_self.gcd(&new_rhs), - } + Self { + poly: self.poly.gcd(&rhs.poly), } } @@ -8622,7 +8316,7 @@ impl PythonMatrix { let mut zero = self.matrix.field.zero(); zero.unify_variables(&mut new_rhs[(0, 0)]); - new_self.field = RationalPolynomialField::new(Z, zero.numerator.get_vars()); + new_self.field = RationalPolynomialField::new(Z); new_rhs.field = new_self.field.clone(); // now update every element @@ -8643,7 +8337,7 @@ impl PythonMatrix { &self, rhs: &PythonRationalPolynomial, ) -> (PythonMatrix, PythonRationalPolynomial) { - if self.matrix.field == RationalPolynomialField::new(Z, rhs.poly.numerator.get_vars()) { + if self.matrix.field == RationalPolynomialField::new(Z) { return (self.clone(), rhs.clone()); } @@ -8653,7 +8347,7 @@ impl PythonMatrix { let mut zero = self.matrix.field.zero(); zero.unify_variables(&mut new_rhs); - new_self.field = RationalPolynomialField::new(Z, zero.numerator.get_vars()); + new_self.field = RationalPolynomialField::new(Z); // now update every element for e in &mut new_self.data { @@ -8679,11 +8373,7 @@ impl PythonMatrix { } Ok(PythonMatrix { - matrix: Matrix::new( - nrows, - ncols, - RationalPolynomialField::new(Z, Arc::new(vec![])), - ), + matrix: Matrix::new(nrows, ncols, RationalPolynomialField::new(Z)), }) } @@ -8697,7 +8387,7 @@ impl PythonMatrix { } Ok(PythonMatrix { - matrix: Matrix::identity(nrows, RationalPolynomialField::new(Z, Arc::new(vec![]))), + matrix: Matrix::identity(nrows, RationalPolynomialField::new(Z)), }) } @@ -8726,7 +8416,7 @@ impl PythonMatrix { } } - let field = RationalPolynomialField::new(Z, first.numerator.get_vars()); + let field = RationalPolynomialField::new(Z); Ok(PythonMatrix { matrix: Matrix::eye(&diag, field), @@ -8758,7 +8448,7 @@ impl PythonMatrix { } } - let field = RationalPolynomialField::new(Z, first.numerator.get_vars()); + let field = RationalPolynomialField::new(Z); Ok(PythonMatrix { matrix: Matrix::new_vec(entries, field), @@ -8792,7 +8482,7 @@ impl PythonMatrix { } } - let field = RationalPolynomialField::new(Z, first.numerator.get_vars()); + let field = RationalPolynomialField::new(Z); Ok(PythonMatrix { matrix: Matrix::from_linear(entries, nrows, ncols, field) @@ -9281,7 +8971,7 @@ impl PythonNumericalIntegrator { match &self.grid { Grid::Continuous(cs) => { let mut a = cs.accumulator.shallow_copy(); - a.update_iter(); + a.update_iter(false); Ok(( a.avg, a.err, @@ -9293,7 +8983,7 @@ impl PythonNumericalIntegrator { } Grid::Discrete(ds) => { let mut a = ds.accumulator.shallow_copy(); - a.update_iter(); + a.update_iter(false); Ok(( a.avg, a.err, diff --git a/src/domains/algebraic_number.rs b/src/domains/algebraic_number.rs index 3790c1fa..3305ad4f 100644 --- a/src/domains/algebraic_number.rs +++ b/src/domains/algebraic_number.rs @@ -41,7 +41,7 @@ where fn to_integer(&self, a: &Self::Element) -> Integer { let mut p = Integer::zero(); for x in a.poly.into_iter() { - p += &(self.poly.field.to_integer(x.coefficient) + p += &(self.poly.ring.to_integer(x.coefficient) * &self.characteristic().pow(x.exponents[0] as u64)); } p @@ -63,7 +63,7 @@ where ::Element: Copy, { AlgebraicExtension::galois_field( - self.poly.field.clone(), + self.poly.ring.clone(), new_pow, self.poly.variables[0].clone(), ) @@ -95,8 +95,8 @@ where 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]); + let (qn, r) = q.quot_rem(&self.poly.ring.size()); + poly.append_monomial(r.to_finite_field(&self.poly.ring), &[pow]); pow += 1; q = qn; } @@ -180,7 +180,7 @@ where { poly.clear(); for (i, c) in coeffs.iter().enumerate() { - poly.append_monomial(poly.field.nth(*c), &[i as u16]); + poly.append_monomial(poly.ring.nth(*c), &[i as u16]); } poly.is_irreducible() @@ -434,7 +434,7 @@ impl Ring for AlgebraicExtension { fn nth(&self, n: u64) -> Self::Element { AlgebraicNumber { - poly: self.poly.constant(self.poly.field.nth(n)), + poly: self.poly.constant(self.poly.ring.nth(n)), } } @@ -459,17 +459,17 @@ impl Ring for AlgebraicExtension { } fn characteristic(&self) -> Integer { - self.poly.field.characteristic() + self.poly.ring.characteristic() } fn size(&self) -> Integer { - self.poly.field.size().pow(self.poly.degree(0) as u64) + self.poly.ring.size().pow(self.poly.degree(0) as u64) } /// 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)) + .map(|_| self.poly.ring.sample(rng, range)) .collect(); let mut poly = self.poly.zero_with_capacity(coeffs.len()); @@ -530,7 +530,7 @@ impl> EuclideanDomain for AlgebraicExtension { let c1 = a.poly.content(); let c2 = b.poly.content(); AlgebraicNumber { - poly: a.poly.constant(a.poly.field.gcd(&c1, &c2)), + poly: a.poly.constant(a.poly.ring.gcd(&c1, &c2)), } } } @@ -573,7 +573,7 @@ impl> AlgebraicExtension { MultivariatePolynomial: Factorize, MultivariatePolynomial>: Factorize, { - assert_eq!(self, &b.field); + assert_eq!(self, &b.ring); let (_, s, g, r) = b.norm_impl(); debug_assert!(r.is_irreducible()); @@ -586,7 +586,7 @@ impl> AlgebraicExtension { 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 y = f.to_element(g2.ring.poly.one().mul_exp(&[1])); let b = f.sub(&y, &f.mul(&a, &f.nth(s as u64))); (f, a, b) @@ -617,12 +617,12 @@ impl, E: Exponent> MultivariatePolynomial, E: Exponent> MultivariatePolynomial FactorizedRationalPolynomialField { poly: &MultivariatePolynomial, ) -> FactorizedRationalPolynomialField { FactorizedRationalPolynomialField { - ring: poly.field.clone(), + ring: poly.ring.clone(), var_map: poly.variables.clone(), _phantom_exp: PhantomData, } @@ -142,8 +142,8 @@ impl FactorizedRationalPolynomial { pub fn is_one(&self) -> bool { self.numerator.is_one() && self.denominators.is_empty() - && self.numerator.field.is_one(&self.numer_coeff) - && self.numerator.field.is_one(&self.denom_coeff) + && self.numerator.ring.is_one(&self.numer_coeff) + && self.numerator.ring.is_one(&self.denom_coeff) } } @@ -158,10 +158,10 @@ impl FromNumeratorAndFactorizedDenominator FactorizedRationalPolynomial { let mut content = num.content(); for (d, _) in &dens { - content = d.field.gcd(&content, &d.content()); + content = d.ring.gcd(&content, &d.content()); } - let (num_int, dens_int) = if num.field.is_one(&content) { + let (num_int, dens_int) = if num.ring.is_one(&content) { ( num.map_coeff(|c| c.numerator(), Z), dens.iter() @@ -170,11 +170,11 @@ impl FromNumeratorAndFactorizedDenominator FromNumeratorAndFactorizedDenominator Self { + if self.get_variables() != other.get_variables() { + let mut a = self.clone(); + let mut b = other.clone(); + a.unify_variables(&mut b); + return a.gcd(&b); + } + let gcd_num = self.numerator.gcd(&other.numerator); let mut disjoint_factors = vec![]; @@ -464,7 +471,7 @@ where } } - let field = &self.numerator.field; + let field = &self.numerator.ring; FactorizedRationalPolynomial { numerator: gcd_num, numer_coeff: field.gcd(&self.numer_coeff, &other.numer_coeff), @@ -588,9 +595,7 @@ where } fn is_one(&self, a: &Self::Element) -> bool { - a.numerator.is_one() - && a.denominators.is_empty() - && a.numerator.field.is_one(&a.denom_coeff) + a.numerator.is_one() && a.denominators.is_empty() && a.numerator.ring.is_one(&a.denom_coeff) } fn one_is_gcd_unit() -> bool { @@ -688,6 +693,13 @@ where return self.clone(); } + if self.get_variables() != other.get_variables() { + let mut a = self.clone(); + let mut b = other.clone(); + a.unify_variables(&mut b); + return &a + &b; + } + let mut den = Vec::with_capacity(self.denominators.len() + other.denominators.len()); let mut num_1 = self.numerator.clone(); let mut num_2 = other.numerator.clone(); @@ -713,7 +725,7 @@ where den.push((d.clone(), *p)); } - let ring = &self.numerator.field; + let ring = &self.numerator.ring; let mut coeff1 = self.numer_coeff.clone(); let mut coeff2 = other.numer_coeff.clone(); let mut new_denom = self.denom_coeff.clone(); @@ -735,8 +747,8 @@ where if num.is_zero() { return FactorizedRationalPolynomial { numerator: num, - numer_coeff: self.numerator.field.zero(), - denom_coeff: self.numerator.field.one(), + numer_coeff: self.numerator.ring.zero(), + denom_coeff: self.numerator.ring.one(), denominators: vec![], }; } @@ -757,11 +769,11 @@ where // make sure the numerator is properly normalized let mut r = - FactorizedRationalPolynomial::from_num_den(num, vec![], &self.numerator.field, false); + FactorizedRationalPolynomial::from_num_den(num, vec![], &self.numerator.ring, false); - let field = &r.numerator.field; + let field = &r.numerator.ring; field.mul_assign(&mut r.numer_coeff, &num_gcd); - let g = r.numerator.field.gcd(&r.numer_coeff, &new_denom); + let g = r.numerator.ring.gcd(&r.numer_coeff, &new_denom); if !field.is_one(&g) { r.numer_coeff = field.quot_rem(&r.numer_coeff, &g).0; new_denom = field.quot_rem(&new_denom, &g).0; @@ -804,7 +816,7 @@ where type Output = Self; fn neg(self) -> Self::Output { FactorizedRationalPolynomial { - numer_coeff: self.numerator.field.neg(&self.numer_coeff), + numer_coeff: self.numerator.ring.neg(&self.numer_coeff), numerator: self.numerator, denom_coeff: self.denom_coeff, denominators: self.denominators, @@ -824,6 +836,13 @@ impl<'a, 'b, R: EuclideanDomain + PolynomialGCD, E: Exponent> return self.clone(); } + if self.get_variables() != other.get_variables() { + let mut a = self.clone(); + let mut b = other.clone(); + a.unify_variables(&mut b); + return &a * &b; + } + let mut reduced_numerator_1 = Cow::Borrowed(&self.numerator); let mut reduced_numerator_2 = Cow::Borrowed(&other.numerator); @@ -870,7 +889,7 @@ impl<'a, 'b, R: EuclideanDomain + PolynomialGCD, E: Exponent> } } - let field = &self.numerator.field; + let field = &self.numerator.ring; let mut constant = field.one(); let mut numer_coeff; @@ -914,6 +933,13 @@ where panic!("Cannot invert 0"); } + if self.get_variables() != other.get_variables() { + let mut a = self.clone(); + let mut b = other.clone(); + a.unify_variables(&mut b); + return &a / &b; + } + let mut reduced_numerator_1 = Cow::Borrowed(&self.numerator); let mut den = Vec::with_capacity(self.denominators.len() + 1); @@ -941,7 +967,7 @@ where let r = FactorizedRationalPolynomial::from_num_den( self.numerator.one(), other.numerator.factor(), - &self.numerator.field, + &self.numerator.ring, false, ); @@ -965,7 +991,7 @@ where } } - let field = &self.numerator.field; + let field = &self.numerator.ring; let denom_coeff = field.mul(&r.denom_coeff, &other.numer_coeff); let mut numer_coeff = other.denom_coeff.clone(); @@ -990,13 +1016,13 @@ where field.mul_assign(&mut numer_coeff, &self.numer_coeff); } let num = reduced_numerator_2 * &reduced_numerator_1; - if !num.field.is_one(&constant) { + if !num.ring.is_one(&constant) { den.push((num.constant(constant), 1)); } // properly normalize the rational polynomial let mut r = - FactorizedRationalPolynomial::from_num_den(num, den, &self.numerator.field, false); + FactorizedRationalPolynomial::from_num_den(num, den, &self.numerator.ring, false); field.mul_assign(&mut r.numer_coeff, &numer_coeff); r } @@ -1033,7 +1059,7 @@ where FactorizedRationalPolynomial::from_num_den( p, vec![], - &self.numerator.field, + &self.numerator.ring, false, ), &exp, @@ -1058,9 +1084,9 @@ where * &FactorizedRationalPolynomial::from_num_den( unfold .numerator - .monomial(self.numerator.field.one(), e.to_vec()), + .monomial(self.numerator.ring.one(), e.to_vec()), vec![], - &self.numerator.field, + &self.numerator.ring, true, )); } @@ -1072,12 +1098,12 @@ where (p.clone(), *pe), (self.numerator.constant(self.denom_coeff.clone()), 1), ], - &self.numerator.field, + &self.numerator.ring, false, ) * &FactorizedRationalPolynomial { numerator: self.numerator.one(), numer_coeff: self.numer_coeff.clone(), - denom_coeff: self.numerator.field.one(), + denom_coeff: self.numerator.ring.one(), denominators: vec![], }); diff --git a/src/domains/rational.rs b/src/domains/rational.rs index 9e68ee2f..759d691e 100644 --- a/src/domains/rational.rs +++ b/src/domains/rational.rs @@ -103,7 +103,7 @@ impl FractionNormalization for Z { impl FractionNormalization for PolynomialRing { fn get_normalization_factor(&self, a: &Self::Element) -> Self::Element { - a.constant(a.field.get_normalization_factor(&a.lcoeff())) + a.constant(a.ring.get_normalization_factor(&a.lcoeff())) } } @@ -1079,9 +1079,9 @@ mod test { f.clone(), ); - let p = PolynomialRing::new_from_poly(&poly2); + let p = PolynomialRing::from_poly(&poly2); let rat = p.to_rational_polynomial(&poly2); - let f = FractionField::new(PolynomialRing::new_from_poly(&rat.numerator)); + let f = FractionField::new(PolynomialRing::from_poly(&rat.numerator)); let b = f.neg(&f.nth(3)); let c = f.add(&rat, &b); diff --git a/src/domains/rational_polynomial.rs b/src/domains/rational_polynomial.rs index 0de54edd..23b2ebad 100644 --- a/src/domains/rational_polynomial.rs +++ b/src/domains/rational_polynomial.rs @@ -27,23 +27,20 @@ use super::{ #[derive(Clone, PartialEq, Eq, Hash, Debug)] pub struct RationalPolynomialField { ring: R, - var_map: Arc>, _phantom_exp: PhantomData, } impl RationalPolynomialField { - pub fn new(coeff_ring: R, var_map: Arc>) -> RationalPolynomialField { + pub fn new(coeff_ring: R) -> RationalPolynomialField { RationalPolynomialField { ring: coeff_ring, - var_map, _phantom_exp: PhantomData, } } - pub fn new_from_poly(poly: &MultivariatePolynomial) -> RationalPolynomialField { + pub fn from_poly(poly: &MultivariatePolynomial) -> RationalPolynomialField { RationalPolynomialField { - ring: poly.field.clone(), - var_map: poly.variables.clone(), + ring: poly.ring.clone(), _phantom_exp: PhantomData, } } @@ -90,7 +87,7 @@ where { fn from(poly: MultivariatePolynomial) -> Self { let d = poly.one(); - let field = poly.field.clone(); + let field = poly.ring.clone(); Self::from_num_den(poly, d, &field, false) } } @@ -110,14 +107,14 @@ where *self = Self::from_num_den( self.numerator.clone(), self.denominator.clone(), - &self.numerator.field, + &self.numerator.ring, false, ); *other = Self::from_num_den( other.numerator.clone(), other.denominator.clone(), - &other.numerator.field, + &other.numerator.ring, false, ); } @@ -177,7 +174,7 @@ impl FromNumeratorAndDenominator field: &IntegerRing, do_gcd: bool, ) -> RationalPolynomial { - let content = num.field.gcd(&num.content(), &den.content()); + let content = num.ring.gcd(&num.content(), &den.content()); let mut num_int = MultivariatePolynomial::new(&Z, None, num.variables); num_int.exponents = num.exponents; @@ -185,7 +182,7 @@ impl FromNumeratorAndDenominator let mut den_int = MultivariatePolynomial::new(&Z, Some(den.nterms()), den.variables); den_int.exponents = den.exponents; - if num.field.is_one(&content) { + if num.ring.is_one(&content) { num_int.coefficients = num .coefficients .into_iter() @@ -200,12 +197,12 @@ impl FromNumeratorAndDenominator num_int.coefficients = num .coefficients .into_iter() - .map(|c| num.field.div(&c, &content).numerator()) + .map(|c| num.ring.div(&c, &content).numerator()) .collect(); den_int.coefficients = den .coefficients .into_iter() - .map(|c| den.field.div(&c, &content).numerator()) + .map(|c| den.ring.div(&c, &content).numerator()) .collect(); } @@ -312,7 +309,7 @@ where panic!("Cannot invert 0"); } - let field = self.numerator.field.clone(); + let field = self.numerator.ring.clone(); Self::from_num_den(self.denominator, self.numerator, &field, false) } @@ -397,7 +394,7 @@ where .append_monomial(e.coefficient.clone(), &e_list_coeff); } else { let mut r = RationalPolynomial::new( - &self.numerator.field.clone(), + &self.numerator.ring.clone(), self.numerator.variables.clone(), ); r.numerator @@ -407,14 +404,14 @@ where } let v = Arc::new(variables.to_vec()); - let field = RationalPolynomialField::new(self.numerator.field.clone(), v.clone()); + let field = RationalPolynomialField::new(self.numerator.ring.clone()); let mut poly = MultivariatePolynomial::new(&field, Some(hm.len()), v); if !ignore_denominator { let denom = RationalPolynomial::from_num_den( self.denominator.one(), self.denominator.clone(), - &self.denominator.field, + &self.denominator.ring, false, ); @@ -433,25 +430,38 @@ where // Convert from a univariate polynomial with rational polynomial coefficients to a rational polynomial. pub fn from_univariate( - f: UnivariatePolynomial>, + mut f: UnivariatePolynomial>, ) -> RationalPolynomial { - let Some(pos) = f - .field - .var_map + if f.is_zero() { + return RationalPolynomial { + numerator: MultivariatePolynomial::new_zero(&f.field.ring), + denominator: MultivariatePolynomial::new_one(&f.field.ring), + }; + } + + let pos = f.coefficients[0] + .get_variables() .iter() .position(|x| x == f.variable.as_ref()) - else { - panic!("Variable not found in the field"); - }; + .unwrap_or_else(|| { + let mut r = RationalPolynomial::new( + &f.coefficients[0].numerator.ring, + Arc::new(vec![f.variable.as_ref().clone()]), + ); + for c in &mut f.coefficients { + c.unify_variables(&mut r); + } + + f.coefficients[0].get_variables().len() - 1 + }); - let mut res = RationalPolynomial::new(&f.field.ring, f.field.var_map.clone()); + let mut res = + RationalPolynomial::new(&f.field.ring, f.coefficients[0].get_variables().clone()); - let mut exp = vec![E::zero(); f.field.var_map.len()]; + let mut exp = vec![E::zero(); f.coefficients[0].get_variables().len()]; exp[pos] = E::one(); - let v: RationalPolynomial = res - .numerator - .monomial(res.numerator.field.one(), exp) - .into(); + let v: RationalPolynomial = + res.numerator.monomial(res.numerator.ring.one(), exp).into(); for (p, c) in f.coefficients.into_iter().enumerate() { res = &res + &(&v.pow(p as u64) * &c); @@ -518,7 +528,7 @@ where } fn zero(&self) -> Self::Element { - let num = MultivariatePolynomial::new(&self.ring, None, self.var_map.clone()); + let num = MultivariatePolynomial::new_zero(&self.ring); RationalPolynomial { denominator: num.one(), numerator: num, @@ -526,7 +536,7 @@ where } fn one(&self) -> Self::Element { - let num = MultivariatePolynomial::new(&self.ring, None, self.var_map.clone()).one(); + let num = MultivariatePolynomial::new_zero(&self.ring).one(); RationalPolynomial { numerator: num.clone(), denominator: num, @@ -645,10 +655,19 @@ where impl<'a, 'b, R: EuclideanDomain + PolynomialGCD + PolynomialGCD, E: Exponent> Add<&'a RationalPolynomial> for &'b RationalPolynomial +where + RationalPolynomial: FromNumeratorAndDenominator, { type Output = RationalPolynomial; fn add(self, other: &'a RationalPolynomial) -> Self::Output { + if self.get_variables() != other.get_variables() { + let mut a = self.clone(); + let mut b = other.clone(); + a.unify_variables(&mut b); + return &a + &b; + } + let denom_gcd = self.denominator.gcd(&other.denominator); let mut a_denom_red = Cow::Borrowed(&self.denominator); @@ -686,7 +705,10 @@ impl<'a, 'b, R: EuclideanDomain + PolynomialGCD + PolynomialGCD, E: Expone } } -impl, E: Exponent> Sub for RationalPolynomial { +impl, E: Exponent> Sub for RationalPolynomial +where + RationalPolynomial: FromNumeratorAndDenominator, +{ type Output = Self; fn sub(self, other: Self) -> Self::Output { @@ -696,11 +718,13 @@ impl, E: Exponent> Sub for RationalPolynom impl<'a, 'b, R: EuclideanDomain + PolynomialGCD, E: Exponent> Sub<&'a RationalPolynomial> for &'b RationalPolynomial +where + RationalPolynomial: FromNumeratorAndDenominator, { type Output = RationalPolynomial; fn sub(self, other: &'a RationalPolynomial) -> Self::Output { - (self.clone()).sub(other.clone()) + self.add(&other.clone().neg()) } } @@ -716,10 +740,19 @@ impl, E: Exponent> Neg for RationalPolynom impl<'a, 'b, R: EuclideanDomain + PolynomialGCD, E: Exponent> Mul<&'a RationalPolynomial> for &'b RationalPolynomial +where + RationalPolynomial: FromNumeratorAndDenominator, { type Output = RationalPolynomial; fn mul(self, other: &'a RationalPolynomial) -> Self::Output { + if self.get_variables() != other.get_variables() { + let mut a = self.clone(); + let mut b = other.clone(); + a.unify_variables(&mut b); + return &a * &b; + } + let gcd1 = self.numerator.gcd(&other.denominator); let gcd2 = self.denominator.gcd(&other.numerator); @@ -781,13 +814,13 @@ where let a = RationalPolynomial::from_num_den( dn, self.denominator.clone(), - &self.numerator.field, + &self.numerator.ring, false, ); let b = RationalPolynomial::from_num_den( &self.numerator * &dd, &self.denominator * &self.denominator, - &self.numerator.field, + &self.numerator.ring, false, ); @@ -806,7 +839,7 @@ where return vec![self.clone()]; } - let rat_field = RationalPolynomialField::new_from_poly(&self.numerator); + let rat_field = RationalPolynomialField::from_poly(&self.numerator); let n = self .numerator .to_univariate(var) @@ -895,7 +928,7 @@ where /// that represents `sum_(r(z) = 0) z*log(a(var, z))` if `r` is dependent on `z`, /// else it is `r*log(a)`. pub fn integrate(&self, var: usize) -> (Vec, Vec<(Self, Self)>) { - let rat_field = RationalPolynomialField::new_from_poly(&self.numerator); + let rat_field = RationalPolynomialField::from_poly(&self.numerator); let n = self .numerator .to_univariate(var) @@ -1005,11 +1038,11 @@ where // create new temporary variable let mut t = MultivariatePolynomial::new( - &self.numerator.field, + &self.numerator.ring, None, Arc::new(vec![Variable::Temporary(0)]), ) - .monomial(self.numerator.field.one(), vec![E::one()]) + .monomial(self.numerator.ring.one(), vec![E::one()]) .into(); let mut w = vec![]; @@ -1024,9 +1057,6 @@ where constant.unify_variables(&mut t); - // adding a variable changes the field - p.field = RationalPolynomialField::new_from_poly(&p.coefficients[0].numerator); - h.field = p.field.clone(); let new_var = p.coefficients[0].numerator.nvars() - 1; let b = h.clone() - p.derivative().mul_coeff(&t); @@ -1072,8 +1102,12 @@ where let ll = RationalPolynomial::from_univariate(r.clone()) .numerator .to_multivariate_polynomial_list(&[var, new_var], true); - let mut bivar_poly = - MultivariatePolynomial::new(&p.field, Some(ll.len()), p.field.var_map.clone()); + + let mut bivar_poly = MultivariatePolynomial::new( + &p.field, + Some(ll.len()), + p.coefficients[0].get_variables().clone(), + ); for (e, p) in ll { bivar_poly.append_monomial(p.into(), &e); } @@ -1117,7 +1151,7 @@ where exp.copy_from_slice(t.exponents); let mm = p.coefficients[0] .numerator - .monomial(self.numerator.field.one(), exp); + .monomial(self.numerator.ring.one(), exp); res = &res + &(t.coefficient * &mm.into()); } @@ -1128,7 +1162,7 @@ where let sol = RationalPolynomial::from_num_den( -a.coefficients[0].clone(), a.coefficients[1].clone(), - &a.coefficients[0].field, + &a.coefficients[0].ring, true, ); @@ -1140,7 +1174,7 @@ where exp.copy_from_slice(t.exponents); let mm = p.coefficients[0] .numerator - .monomial(self.numerator.field.one(), exp); + .monomial(self.numerator.ring.one(), exp); res = &res + &(t.coefficient * &mm.into()); } @@ -1172,7 +1206,7 @@ mod test { #[test] fn field() { - let field = RationalPolynomialField::<_, u8>::new(Z, Arc::new(vec![])); + let field = RationalPolynomialField::<_, u8>::new(Z); let one = field.one(); let t = format!("{}", field.printer(&one)); assert_eq!(t, "1"); diff --git a/src/numerical_integration.rs b/src/numerical_integration.rs index de5ea982..1d17be07 100644 --- a/src/numerical_integration.rs +++ b/src/numerical_integration.rs @@ -8,13 +8,13 @@ use crate::domains::float::{ConstructibleFloat, Real, RealNumberLike}; /// the error and the chi-squared of samples added over multiple /// iterations. /// -/// Samples can be added using [`StatisticsAccumulator::add_samples()`]. When an iteration of -/// samples is finished, call [`StatisticsAccumulator::update_iter()`], which +/// Samples can be added using [`Self::add_sample()`]. When an iteration of +/// samples is finished, call [`Self::update_iter()`], which /// updates the average, error and chi-squared over all iterations with the average /// and error of the current iteration in a weighted fashion. /// -/// This accumulator can be merged with another accumulator using [`StatisticsAccumulator::merge()`] or -/// [`StatisticsAccumulator::merge_samples_no_reset()`]. This is useful when +/// This accumulator can be merged with another accumulator using [`Self::merge_samples()`] or +/// [`Self::merge_samples_no_reset()`]. This is useful when /// samples are collected in multiple threads. /// /// The accumulator also stores which samples yielded the highest weight thus far. @@ -24,6 +24,8 @@ pub struct StatisticsAccumulator Statisti StatisticsAccumulator { sum: T::new_zero(), sum_sq: T::new_zero(), + total_sum: T::new_zero(), + total_sum_sq: T::new_zero(), weight_sum: T::new_zero(), avg_sum: T::new_zero(), avg: T::new_zero(), @@ -77,6 +81,8 @@ impl Statisti StatisticsAccumulator { sum: self.sum, sum_sq: self.sum_sq, + total_sum: self.total_sum, + total_sum_sq: self.total_sum_sq, weight_sum: self.weight_sum, avg_sum: self.avg_sum, avg: self.avg, @@ -152,9 +158,12 @@ impl Statisti } } - /// Process the samples added with `[`Self::add_sample()`]` and + /// Process the samples added with [`Self::add_sample()`] and /// compute a new average, error, and chi-squared. - pub fn update_iter(&mut self) -> bool { + /// + /// When `weighted_average=True`, a weighted average and error is computed using + /// the iteration variances as a weight. + pub fn update_iter(&mut self, weighted_average: bool) -> bool { // TODO: we could be throwing away events that are very rare if self.new_samples < 2 { self.cur_iter += 1; @@ -164,9 +173,11 @@ impl Statisti self.processed_samples += self.new_samples; self.num_zero_evaluations += self.new_zero_evaluations; let n = T::new_from_usize(self.new_samples); - self.sum /= &n; - self.sum_sq /= n * n; - let mut w = (self.sum_sq * n).sqrt(); + self.total_sum += self.sum; + self.total_sum_sq += self.sum_sq; + self.sum /= n; + self.sum_sq /= n; + let mut w = self.sum_sq.sqrt(); w = ((w + self.sum) * (w - self.sum)) / (n - T::new_one()); if w == T::new_zero() { @@ -179,9 +190,17 @@ impl Statisti self.weight_sum += w; self.avg_sum += w * self.sum; - let sigma_sq = self.weight_sum.inv(); - self.avg = sigma_sq * self.avg_sum; - self.err = sigma_sq.sqrt(); + + if weighted_average { + let sigma_sq = self.weight_sum.inv(); + self.avg = sigma_sq * self.avg_sum; + self.err = sigma_sq.sqrt(); + } else { + let n = T::new_from_usize(self.processed_samples); + self.avg = self.total_sum / n; + self.err = ((self.total_sum_sq / n - self.avg * self.avg) / (n - T::new_one())).sqrt(); + } + if self.cur_iter == 0 { self.guess = self.sum; } @@ -200,6 +219,23 @@ impl Statisti true } + /// Get an estimate for the average, error and chi-squared, as if the current iteration + /// has ended without adding more samples. + pub fn get_live_estimate(&self) -> (T, T, T) { + let mut a = self.shallow_copy(); + a.update_iter(false); + (a.avg, a.err, a.chi_sq) + } + + /// Format the live `mean ± sdev` as `mean(sdev)` in a human-readable way with the correct number of digits. + /// + /// Based on the Python package [gvar](https://github.com/gplepage/gvar) by Peter Lepage. + pub fn format_live_uncertainty(&self) -> String { + let mut a = self.shallow_copy(); + a.update_iter(false); + Self::format_uncertainty_impl(a.avg.to_f64(), a.err.to_f64()) + } + /// Format `mean ± sdev` as `mean(sdev)` in a human-readable way with the correct number of digits. /// /// Based on the Python package [gvar](https://github.com/gplepage/gvar) by Peter Lepage. @@ -536,7 +572,7 @@ impl Discrete } let acc = &mut bin.accumulator; - acc.update_iter(); + acc.update_iter(false); if acc.processed_samples > 1 { err_sum += acc.err * T::new_from_usize(acc.processed_samples - 1).sqrt(); @@ -587,7 +623,7 @@ impl Discrete bin.pdf /= sum; } - self.accumulator.update_iter(); + self.accumulator.update_iter(false); } /// Sample a point form this grid, writing the result in `sample`. @@ -772,7 +808,7 @@ impl Continuo d.update(learning_rate); } - self.accumulator.update_iter(); + self.accumulator.update_iter(false); } /// Returns `Ok` when this grid can be merged with another grid, @@ -1133,7 +1169,7 @@ mod test { grid.update(1.5, 1.5); } - assert_eq!(grid.accumulator.avg, 0.9713543844460519); - assert_eq!(grid.accumulator.err, 0.0009026050146732183) + assert_eq!(grid.accumulator.avg, 0.9718412953459551); + assert_eq!(grid.accumulator.err, 0.0009349254838085983) } } diff --git a/src/poly.rs b/src/poly.rs index 692cb8b6..d16411cd 100644 --- a/src/poly.rs +++ b/src/poly.rs @@ -653,7 +653,7 @@ impl<'a> AtomView<'a> { poly: &mut MultivariatePolynomial, field: &R, ) { - let mut coefficient = poly.field.one(); + let mut coefficient = poly.ring.one(); let mut exponents = smallvec![E::zero(); vars.len()]; match term { @@ -1221,7 +1221,7 @@ impl MultivariatePolynomial { } } - f(&self.field, monomial.coefficient, &mut coeff); + f(&self.ring, monomial.coefficient, &mut coeff); mul.extend(coeff.as_view()); add.extend(mul_h.as_view()); } @@ -1355,7 +1355,7 @@ impl Token { poly: &mut MultivariatePolynomial, field: &R, ) -> Result<(), Cow<'static, str>> { - let mut coefficient = poly.field.one(); + let mut coefficient = poly.ring.one(); let mut exponents = smallvec![E::zero(); var_name_map.len()]; match term { diff --git a/src/poly/evaluate.rs b/src/poly/evaluate.rs index 8d13af3f..b95a1ac3 100644 --- a/src/poly/evaluate.rs +++ b/src/poly/evaluate.rs @@ -351,7 +351,7 @@ impl MultivariatePolynomial { debug_assert!(indices.len() <= index_start + 1); let num = if indices.len() == index_start { - self.field.zero() + self.ring.zero() } else { self.coefficients[indices[index_start]].clone() }; diff --git a/src/poly/factor.rs b/src/poly/factor.rs index feeca181..c8531dc8 100644 --- a/src/poly/factor.rs +++ b/src/poly/factor.rs @@ -343,7 +343,7 @@ impl Factorize for MultivariatePolynomial Factorize for MultivariatePolynomial Factorize for MultivariatePolynomial Factorize let mut factors = vec![]; - if !self.field.is_one(&c) { + if !self.ring.is_one(&c) { factors.push((self.constant(c), 1)); } @@ -455,15 +455,15 @@ impl Factorize return vec![(f.clone(), 1)]; } - let mut g_f = g.to_number_field(&self.field); + let mut g_f = g.to_number_field(&self.ring); for (f, b) in factors { debug!("Rational factor {}", f); let alpha_poly = g.variable(&self.get_vars_ref()[v]).unwrap() - + g.variable(&self.field.poly().variables[0]).unwrap() + + g.variable(&self.ring.poly().variables[0]).unwrap() * &g.constant((s as u64).into()); - let f = f.to_number_field(&self.field); + let f = f.to_number_field(&self.ring); let gcd = f.gcd(&g_f); @@ -471,7 +471,7 @@ impl Factorize let g = MultivariatePolynomial::from_number_field(&gcd) .replace_with_poly(v, &alpha_poly); - full_factors.push((g.to_number_field(&self.field), b * p)); + full_factors.push((g.to_number_field(&self.ring), b * p)); } } @@ -506,7 +506,7 @@ where factors.append(&mut nf); } - if factors.is_empty() || !self.field.is_one(&c) { + if factors.is_empty() || !self.ring.is_one(&c) { factors.push((self.constant(c), 1)) } @@ -680,7 +680,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.characteristic().to_u64() as usize; + let p = self.ring.characteristic().to_u64() as usize; let mut b = f.clone(); for es in b.exponents_iter_mut() { for e in es { @@ -738,7 +738,7 @@ where let mut factors = vec![]; let mut i = 1; - while !w.is_constant() && i < self.field.characteristic().to_u64() as usize { + while !w.is_constant() && i < self.ring.characteristic().to_u64() as usize { let z = v - w.derivative(var); let g = w.gcd(&z); w = w / &g; @@ -762,7 +762,7 @@ where let mut e = self.last_exponents().to_vec(); e[var] = E::one(); - let x = self.monomial(self.field.one(), e); + let x = self.monomial(self.ring.one(), e); let mut factors = vec![]; let mut h = x.clone(); @@ -771,7 +771,7 @@ where while !f.is_one() { i += 1; - h = h.exp_mod_univariate(self.field.size(), &mut f); + h = h.exp_mod_univariate(self.ring.size(), &mut f); let mut g = f.gcd(&(&h - &x)); @@ -816,7 +816,7 @@ where let mut exp = vec![E::zero(); self.nvars()]; let mut try_counter = 0; - let characteristic = self.field.characteristic(); + let characteristic = self.ring.characteristic(); let factor = loop { // generate a random non-constant polynomial @@ -824,14 +824,14 @@ where if d == 1 && (characteristic.is_zero() || try_counter < characteristic) { exp[var] = E::zero(); - random_poly.append_monomial(self.field.nth(try_counter), &exp); + random_poly.append_monomial(self.ring.nth(try_counter), &exp); exp[var] = E::one(); - random_poly.append_monomial(self.field.one(), &exp); + random_poly.append_monomial(self.ring.one(), &exp); try_counter += 1; } else { for i in 0..2 * d { let r = self - .field + .ring .sample(&mut rng, (0, characteristic.to_i64().unwrap_or(i64::MAX))); if !F::is_zero(&r) { exp[var] = E::from_u32(i as u32); @@ -850,8 +850,8 @@ where break g; } - let b = if self.field.characteristic() == 2 { - let max = self.field.get_extension_degree() as usize * d; + let b = if self.ring.characteristic() == 2 { + let max = self.ring.get_extension_degree() as usize * d; let mut b = random_poly.clone(); let mut vcur = b.clone(); @@ -864,7 +864,7 @@ where b } else { // TODO: use Frobenius map and modular composition to prevent computing large exponent poly^(p^d) - let p = self.field.size(); + let p = self.ring.size(); random_poly .exp_mod_univariate(&(&p.pow(d as u64) - &1i64.into()) / &2i64.into(), &mut s) - self.one() @@ -910,7 +910,7 @@ where let y_poly = self.to_univariate_polynomial_list(interpolation_var); // add the leading coefficient as a first factor - let mut factors = vec![lcoeff.replace(interpolation_var, &self.field.zero())]; + let mut factors = vec![lcoeff.replace(interpolation_var, &self.ring.zero())]; factors.extend_from_slice(univariate_factors); // extract coefficients in y @@ -1009,29 +1009,29 @@ where return factors; } - let mut sample_point = self.field.zero(); + let mut sample_point = self.ring.zero(); let mut uni_f = self.replace(interpolation_var, &sample_point); let mut i = 0; let mut rng = thread_rng(); loop { - if self.field.size() == i { + if self.ring.size() == i { let field = self - .field - .upgrade(self.field.get_extension_degree().to_u64() as usize + 1); + .ring + .upgrade(self.ring.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 s_l = self.map_coeff(|c| self.ring.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())) + .map(|f| f.map_coeff(|c| self.ring.downgrade_element(c), self.ring.clone())) .collect(); } @@ -1042,7 +1042,7 @@ where } // TODO: sample simple points first - sample_point = self.field.sample(&mut rng, (0, i)); + sample_point = self.ring.sample(&mut rng, (0, i)); uni_f = self.replace(interpolation_var, &sample_point); i += 1; } @@ -1121,7 +1121,7 @@ where 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)); + *x = x.shift_var_cached(interpolation_var, &self.ring.neg(&sample_point)); } } @@ -1131,7 +1131,7 @@ where /// Reconstruct the leading coefficient using a Pade approximation with numerator degree `deg_n` and /// denominator degree `deg_d`. The resulting denominator should be a factor of the leading coefficient. fn lcoeff_reconstruct(coeffs: &[Self], deg_n: u32, deg_d: u32) -> Self { - let mut lcoeff = coeffs[0].constant(coeffs[0].field.one()); + let mut lcoeff = coeffs[0].constant(coeffs[0].ring.one()); for x in coeffs { let d = x.rational_approximant_univariate(deg_n, deg_d).unwrap().1; if !d.is_one() { @@ -1421,7 +1421,7 @@ where let univariate_deltas = Self::diophantine_univariate( &mut univariate_factors, - &factors[0].constant(factors[0].field.one()), + &factors[0].constant(factors[0].ring.one()), ); (univariate_factors, univariate_deltas) @@ -1467,7 +1467,7 @@ where // select a suitable evaluation point let mut sample_points: Vec<_> = - order[1..].iter().map(|i| (*i, self.field.zero())).collect(); + order[1..].iter().map(|i| (*i, self.ring.zero())).collect(); let mut uni_f; let mut biv_f; let mut rng = thread_rng(); @@ -1480,18 +1480,18 @@ where 'new_sample: loop { sample_fail += &1.into(); - if &sample_fail * &2.into() > self.field.size() { + if &sample_fail * &2.into() > self.ring.size() { // the field is too small, upgrade let field = self - .field - .upgrade(self.field.get_extension_degree().to_u64() as usize + 1); + .ring + .upgrade(self.ring.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 s_l = self.map_coeff(|c| self.ring.upgrade_element(c, &field), field.clone()); let facs = s_l.multivariate_factorization( order, @@ -1501,12 +1501,12 @@ where return facs .into_iter() - .map(|f| f.map_coeff(|c| self.field.downgrade_element(c), self.field.clone())) + .map(|f| f.map_coeff(|c| self.ring.downgrade_element(c), self.ring.clone())) .collect(); } for s in &mut sample_points { - s.1 = self.field.nth(rng.gen_range(0..=coefficient_upper_bound)); + s.1 = self.ring.nth(rng.gen_range(0..=coefficient_upper_bound)); } biv_f = self.clone(); @@ -1566,7 +1566,7 @@ where } for (v, s) in &sample_points { - debug!("Sample point {}={}", v, self.field.printer(s)); + debug!("Sample point {}={}", v, self.ring.printer(s)); } let bivariate_factors = biv_f.bivariate_factorization(order[0], order[1]); @@ -1692,8 +1692,8 @@ impl MultivariatePolynomial { let shift = &sample_points.iter().find(|s| s.0 == *r.0).unwrap().1; exp[*r.0] = E::one(); let var_pow = error - .monomial(error.field.one(), exp.clone()) - .shift_var(*r.0, &error.field.neg(shift)) + .monomial(error.ring.one(), exp.clone()) + .shift_var(*r.0, &error.ring.neg(shift)) .pow(r.1 + 1); exp[*r.0] = E::zero(); mod_vars.push(var_pow); @@ -1701,8 +1701,8 @@ impl MultivariatePolynomial { exp[last_var] = E::one(); let var_pow = error - .monomial(error.field.one(), exp) - .shift_var(last_var, &error.field.neg(shift)); + .monomial(error.ring.one(), exp) + .shift_var(last_var, &error.ring.neg(shift)); let mut cur_exponent; let mut next_exponent = var_pow.clone(); @@ -1856,7 +1856,7 @@ impl MultivariatePolynomial { ); for f in &mut reconstructed_factors { - *f = f.shift_var(order[v], &self.field.neg(shift)); + *f = f.shift_var(order[v], &self.ring.neg(shift)); } for f in &reconstructed_factors { @@ -1909,7 +1909,7 @@ impl MultivariatePolynomial { let prod_mod = prods .iter() - .map(|f| f.replace(last_var, &self.field.zero())) + .map(|f| f.replace(last_var, &self.ring.zero())) .collect::>(); debug!("in shift {}", self); @@ -2003,9 +2003,9 @@ impl MultivariatePolynomial { { let lcoeff = self.lcoeff(); // lcoeff % p != 0 let mut gamma = gamma.unwrap_or(lcoeff.clone()); - let lcoeff_p = lcoeff.to_finite_field(&u.field); - let gamma_p = gamma.to_finite_field(&u.field); - let field = u.field.clone(); + let lcoeff_p = lcoeff.to_finite_field(&u.ring); + let gamma_p = gamma.to_finite_field(&u.ring); + let field = u.ring.clone(); let p = Integer::from(field.get_prime().to_u64()); let a = self.clone().mul_coeff(gamma.clone()); @@ -2253,7 +2253,7 @@ impl MultivariatePolynomial { let finite_field = Zp::new(p); // add the leading coefficient as a first factor - let mut factors = vec![lcoeff.replace(interpolation_var, &poly.field.zero())]; + let mut factors = vec![lcoeff.replace(interpolation_var, &poly.ring.zero())]; for f in univariate_factors { factors.push(f.clone()); @@ -2261,7 +2261,7 @@ impl MultivariatePolynomial { let delta = Self::lift_diophantine_univariate( &mut factors, - &poly.constant(poly.field.one()), + &poly.constant(poly.ring.one()), finite_field.get_prime(), k, ); @@ -2511,7 +2511,7 @@ impl MultivariatePolynomial { if !sample_point.is_zero() { for x in &mut rec_factors { // shift the polynomial to y - sample - *x = x.shift_var(interpolation_var, &self.field.neg(&sample_point)); + *x = x.shift_var(interpolation_var, &self.ring.neg(&sample_point)); } } @@ -2543,8 +2543,8 @@ impl MultivariatePolynomial { .iter() .map(|s| { s.map_coeff( - |c| field.to_symmetric_integer(c).to_finite_field(&rhs.field), - rhs.field.clone(), + |c| field.to_symmetric_integer(c).to_finite_field(&rhs.ring), + rhs.ring.clone(), ) }) .collect(); @@ -2553,7 +2553,7 @@ impl MultivariatePolynomial { return deltas; } - let mut tot = rhs.constant(rhs.field.one()); + let mut tot = rhs.constant(rhs.ring.one()); for f in factors.iter() { tot = &tot * f; } @@ -2579,8 +2579,8 @@ impl MultivariatePolynomial { *d = &*d + &new_delta.map_coeff( - |c| (&field.to_symmetric_integer(c) * &m).to_finite_field(&rhs.field), - rhs.field.clone(), + |c| (&field.to_symmetric_integer(c) * &m).to_finite_field(&rhs.ring), + rhs.ring.clone(), ); } @@ -3359,7 +3359,7 @@ impl MultivariatePolynomial, E, LexOrder> { let univariate_deltas = MultivariatePolynomial::lift_diophantine_univariate( &mut univariate_factors, - &factors[0].constant(factors[0].field.one()), + &factors[0].constant(factors[0].ring.one()), p, k, ); diff --git a/src/poly/gcd.rs b/src/poly/gcd.rs index 0bec5dcd..838ddca7 100755 --- a/src/poly/gcd.rs +++ b/src/poly/gcd.rs @@ -130,7 +130,7 @@ impl MultivariatePolynomial { r: &[(usize, R::Element)], cache: &mut [Vec], ) -> Vec { - let mut eval = vec![self.field.one(); self.nterms()]; + let mut eval = vec![self.ring.one(); self.nterms()]; for (c, t) in eval.iter_mut().zip(self) { // evaluate each exponent for (n, v) in r { @@ -138,12 +138,12 @@ impl MultivariatePolynomial { if exp > 0 { if exp < cache[*n].len() { if R::is_zero(&cache[*n][exp]) { - cache[*n][exp] = self.field.pow(v, exp as u64); + cache[*n][exp] = self.ring.pow(v, exp as u64); } - self.field.mul_assign(c, &cache[*n][exp]); + self.ring.mul_assign(c, &cache[*n][exp]); } else { - self.field.mul_assign(c, &self.field.pow(v, exp as u64)); + self.ring.mul_assign(c, &self.ring.pow(v, exp as u64)); } } } @@ -160,7 +160,7 @@ impl MultivariatePolynomial { out: &mut MultivariatePolynomial, ) { out.clear(); - let mut c = self.field.zero(); + let mut c = self.ring.zero(); let mut new_exp = vec![E::zero(); self.nvars()]; for (aa, e) in self.into_iter().zip(exp_evals) { if aa.exponents[main_var] != new_exp[main_var] { @@ -168,13 +168,13 @@ impl MultivariatePolynomial { out.coefficients.push(c); out.exponents.extend_from_slice(&new_exp); - c = self.field.zero(); + c = self.ring.zero(); } new_exp[main_var] = aa.exponents[main_var]; } - self.field.add_mul_assign(&mut c, aa.coefficient, e); + self.ring.add_mul_assign(&mut c, aa.coefficient, e); } if !R::is_zero(&c) { @@ -212,7 +212,7 @@ impl MultivariatePolynomial { // normalize the gcd let l = d.coefficients.last().unwrap().clone(); for x in &mut d.coefficients { - self.field.div_assign(x, &l); + self.ring.div_assign(x, &l); } d @@ -234,19 +234,19 @@ impl MultivariatePolynomial { if exp > 0 { if exp < cache[*n].len() { if F::is_zero(&cache[*n][exp]) { - cache[*n][exp] = self.field.pow(vv, exp as u64); + cache[*n][exp] = self.ring.pow(vv, exp as u64); } - self.field.mul_assign(&mut c, &cache[*n][exp]); + self.ring.mul_assign(&mut c, &cache[*n][exp]); } else { - self.field - .mul_assign(&mut c, &self.field.pow(vv, exp as u64)); + self.ring + .mul_assign(&mut c, &self.ring.pow(vv, exp as u64)); } } } tm.entry(mv.exponents[v]) - .and_modify(|e| self.field.add_assign(e, &c)) + .and_modify(|e| self.ring.add_assign(e, &c)) .or_insert(c); } @@ -275,7 +275,7 @@ impl MultivariatePolynomial { let mut cache = (0..ap.nvars()) .map(|i| { vec![ - ap.field.zero(); + ap.ring.zero(); min( max(ap.degree(i), bp.degree(i)).to_u32() as usize + 1, POW_CACHE_SIZE @@ -295,13 +295,13 @@ impl MultivariatePolynomial { let (_, a1, b1) = loop { for v in &mut cache { for vi in v { - *vi = ap.field.zero(); + *vi = ap.ring.zero(); } } let r: Vec<_> = vars .iter() - .map(|i| (*i, ap.field.sample(&mut rng, (1, MAX_RNG_PREFACTOR as i64)))) + .map(|i| (*i, ap.ring.sample(&mut rng, (1, MAX_RNG_PREFACTOR as i64)))) .collect(); let a1 = ap.sample_polynomial(var, &r, &mut cache, &mut tm); @@ -311,7 +311,7 @@ impl MultivariatePolynomial { break (r, a1, b1); } - if !ap.field.size().is_zero() && fail_count * 2 > ap.field.size() { + if !ap.ring.size().is_zero() && fail_count * 2 > ap.ring.size() { debug!("Field is too small to find a good sample point"); // TODO: upgrade to larger field? return ap.degree(var).min(bp.degree(var)); @@ -340,7 +340,7 @@ impl MultivariatePolynomial { // solve the transposed Vandermonde system for (((c, ex), sample), rhs) in shape.iter().zip(&row_sample_values).zip(&samples) { if c.nterms() == 1 { - let coeff = self.field.div(&rhs[0], &sample[0]); + let coeff = self.ring.div(&rhs[0], &sample[0]); let mut ee: SmallVec<[E; INLINED_EXPONENTS]> = c.exponents(0).into(); ee[main_var] = *ex; gp.append_monomial(coeff, &ee); @@ -348,43 +348,43 @@ impl MultivariatePolynomial { } // construct the master polynomial (1-s1)*(1-s2)*... efficiently - let mut master = vec![self.field.zero(); sample.len() + 1]; - master[0] = self.field.one(); + let mut master = vec![self.ring.zero(); sample.len() + 1]; + master[0] = self.ring.one(); for (i, x) in sample.iter().take(c.nterms()).enumerate() { let first = &mut master[0]; let mut old_last = first.clone(); - self.field.mul_assign(first, &self.field.neg(x)); + self.ring.mul_assign(first, &self.ring.neg(x)); for m in &mut master[1..=i] { let ov = m.clone(); - self.field.mul_assign(m, &self.field.neg(x)); - self.field.add_assign(m, &old_last); + self.ring.mul_assign(m, &self.ring.neg(x)); + self.ring.add_assign(m, &old_last); old_last = ov; } - master[i + 1] = self.field.one(); + master[i + 1] = self.ring.one(); } for (i, s) in sample.iter().take(c.nterms()).enumerate() { - let mut norm = self.field.one(); + let mut norm = self.ring.one(); // sample master/(1-s_i) by using the factorized form for (j, l) in sample.iter().enumerate() { if j != i { - self.field.mul_assign(&mut norm, &self.field.sub(s, l)) + self.ring.mul_assign(&mut norm, &self.ring.sub(s, l)) } } // divide out 1-s_i - let mut coeff = self.field.zero(); - let mut last_q = self.field.zero(); + let mut coeff = self.ring.zero(); + let mut last_q = self.ring.zero(); for (m, rhs) in master.iter().skip(1).zip(rhs).rev() { - last_q = self.field.add(m, &self.field.mul(s, &last_q)); - self.field.add_mul_assign(&mut coeff, &last_q, rhs); + last_q = self.ring.add(m, &self.ring.mul(s, &last_q)); + self.ring.add_mul_assign(&mut coeff, &last_q, rhs); } - self.field.div_assign(&mut coeff, &norm); + self.ring.div_assign(&mut coeff, &norm); // divide by the Vandermonde row since the Vandermonde matrices should start with a 1 - self.field.div_assign(&mut coeff, s); + self.ring.div_assign(&mut coeff, s); let mut ee: SmallVec<[E; INLINED_EXPONENTS]> = c.exponents(i).into(); ee[main_var] = *ex; @@ -403,16 +403,16 @@ impl MultivariatePolynomial { u: &[MultivariatePolynomial], x: usize, // the variable index to extend the polynomial by ) -> MultivariatePolynomial { - let field = &u[0].field; + let field = &u[0].ring; // compute inverses let mut gammas = Vec::with_capacity(a.len()); for k in 1..a.len() { let mut pr = field.sub(&a[k], &a[0]); for i in 1..k { - u[0].field.mul_assign(&mut pr, &field.sub(&a[k], &a[i])); + u[0].ring.mul_assign(&mut pr, &field.sub(&a[k], &a[i])); } - gammas.push(u[0].field.inv(&pr)); + gammas.push(u[0].ring.inv(&pr)); } // compute Newton coefficients @@ -481,7 +481,7 @@ impl MultivariatePolynomial { let (_, d) = &shape[single_scale]; for t in &g { if t.exponents[main_var] == *d { - let scale_factor = a.field.neg(&a.field.inv(t.coefficient)); // TODO: why -1? + let scale_factor = a.ring.neg(&a.ring.inv(t.coefficient)); // TODO: why -1? return Ok(g.mul_coeff(scale_factor)); } } @@ -499,7 +499,7 @@ impl MultivariatePolynomial { let mut cache = (0..a.nvars()) .map(|i| { vec![ - a.field.zero(); + a.ring.zero(); min( max(a.degree(i), b.degree(i)).to_u32() as usize + 1, POW_CACHE_SIZE @@ -512,13 +512,13 @@ impl MultivariatePolynomial { let (row_sample_values, samples) = 'find_root_sample: loop { for v in &mut cache { for vi in v { - *vi = a.field.zero(); + *vi = a.ring.zero(); } } let r_orig: SmallVec<[_; INLINED_EXPONENTS]> = vars .iter() - .map(|i| (*i, a.field.sample(&mut rng, (1, MAX_RNG_PREFACTOR as i64)))) + .map(|i| (*i, a.ring.sample(&mut rng, (1, MAX_RNG_PREFACTOR as i64)))) .collect(); let mut row_sample_values = Vec::with_capacity(shape.len()); // coefficients for the linear system @@ -530,18 +530,18 @@ impl MultivariatePolynomial { for t in c { // evaluate each exponent - let mut c = a.field.one(); + let mut c = a.ring.one(); for (n, v) in &r_orig { let exp = t.exponents[*n].to_u32() as usize; if exp > 0 { if exp < cache[*n].len() { if F::is_zero(&cache[*n][exp]) { - cache[*n][exp] = a.field.pow(v, exp as u64); + cache[*n][exp] = a.ring.pow(v, exp as u64); } - a.field.mul_assign(&mut c, &cache[*n][exp]); + a.ring.mul_assign(&mut c, &cache[*n][exp]); } else { - a.field.mul_assign(&mut c, &a.field.pow(v, exp as u64)); + a.ring.mul_assign(&mut c, &a.ring.pow(v, exp as u64)); } } } @@ -573,14 +573,14 @@ impl MultivariatePolynomial { // sample at r^i if sample_index > 0 { for (c, rr) in r.iter_mut().zip(&r_orig) { - *c = (c.0, a.field.mul(&c.1, &rr.1)); + *c = (c.0, a.ring.mul(&c.1, &rr.1)); } for (c, e) in a_current.to_mut().iter_mut().zip(&a_eval) { - a.field.mul_assign(c, e); + a.ring.mul_assign(c, e); } for (c, e) in b_current.to_mut().iter_mut().zip(&b_eval) { - b.field.mul_assign(c, e); + b.ring.mul_assign(c, e); } } @@ -613,7 +613,7 @@ impl MultivariatePolynomial { // p is likely unlucky debug!( "Bad current image: gcd({},{}) mod {} under {:?} = {}", - a, b, a.field, r, g + a, b, a.ring, r, g ); return Err(GCDError::BadCurrentImage); } @@ -622,21 +622,21 @@ impl MultivariatePolynomial { } // construct the scaling coefficient - let mut scale_factor = a.field.one(); - let mut coeff = a.field.one(); + let mut scale_factor = a.ring.one(); + let mut coeff = a.ring.one(); let (c, d) = &shape[single_scale]; for (n, v) in r.iter() { // TODO: can be taken from row? - a.field.mul_assign( + a.ring.mul_assign( &mut coeff, - &a.field.pow(v, c.exponents(0)[*n].to_u32() as u64), + &a.ring.pow(v, c.exponents(0)[*n].to_u32() as u64), ); } let mut found = false; for t in &g { if t.exponents[main_var] == *d { - scale_factor = g.field.div(&coeff, t.coefficient); + scale_factor = g.ring.div(&coeff, t.coefficient); found = true; break; } @@ -667,17 +667,17 @@ impl MultivariatePolynomial { // find the associated term in the sample, trying the usual place first if i < g.nterms() && g.exponents(i)[main_var] == *exp { - rhs.push(a.field.neg(&a.field.mul(&g.coefficients[i], &scale_factor))); + rhs.push(a.ring.neg(&a.ring.mul(&g.coefficients[i], &scale_factor))); } else { // find the matching term if it exists for m in g.into_iter() { if m.exponents[main_var] == *exp { - rhs.push(a.field.neg(&a.field.mul(m.coefficient, &scale_factor))); + rhs.push(a.ring.neg(&a.ring.mul(m.coefficient, &scale_factor))); continue 'rhs; } } - rhs.push(a.field.zero()); + rhs.push(a.ring.zero()); } } } @@ -711,7 +711,7 @@ impl MultivariatePolynomial { let mut cache = (0..a.nvars()) .map(|i| { vec![ - a.field.zero(); + a.ring.zero(); min( max(a.degree(i), b.degree(i)).to_u32() as usize + 1, POW_CACHE_SIZE @@ -732,13 +732,13 @@ impl MultivariatePolynomial { let (row_sample_values, samples) = 'find_root_sample: loop { for v in &mut cache { for vi in v { - *vi = a.field.zero(); + *vi = a.ring.zero(); } } let r_orig: SmallVec<[_; INLINED_EXPONENTS]> = vars .iter() - .map(|i| (*i, a.field.sample(&mut rng, (1, MAX_RNG_PREFACTOR as i64)))) + .map(|i| (*i, a.ring.sample(&mut rng, (1, MAX_RNG_PREFACTOR as i64)))) .collect(); let mut row_sample_values = Vec::with_capacity(shape.len()); // coefficients for the linear system @@ -750,18 +750,18 @@ impl MultivariatePolynomial { for t in c { // evaluate each exponent - let mut c = a.field.one(); + let mut c = a.ring.one(); for (n, v) in &r_orig { let exp = t.exponents[*n].to_u32() as usize; if exp > 0 { if exp < cache[*n].len() { if F::is_zero(&cache[*n][exp]) { - cache[*n][exp] = a.field.pow(v, exp as u64); + cache[*n][exp] = a.ring.pow(v, exp as u64); } - a.field.mul_assign(&mut c, &cache[*n][exp]); + a.ring.mul_assign(&mut c, &cache[*n][exp]); } else { - a.field.mul_assign(&mut c, &a.field.pow(v, exp as u64)); + a.ring.mul_assign(&mut c, &a.ring.pow(v, exp as u64)); } } } @@ -800,14 +800,14 @@ impl MultivariatePolynomial { // sample at r^i if sample_index > 0 { for (c, rr) in r.iter_mut().zip(&r_orig) { - *c = (c.0, a.field.mul(&c.1, &rr.1)); + *c = (c.0, a.ring.mul(&c.1, &rr.1)); } for (c, e) in a_current.to_mut().iter_mut().zip(&a_eval) { - a.field.mul_assign(c, e); + a.ring.mul_assign(c, e); } for (c, e) in b_current.to_mut().iter_mut().zip(&b_eval) { - b.field.mul_assign(c, e); + b.ring.mul_assign(c, e); } } @@ -840,7 +840,7 @@ impl MultivariatePolynomial { // p is likely unlucky debug!( "Bad current image: gcd({},{}) mod {} under {:?} = {}", - a, b, a.field, r, g + a, b, a.ring, r, g ); return Err(GCDError::BadCurrentImage); } @@ -862,7 +862,7 @@ impl MultivariatePolynomial { let mut found = false; for t in &g { if t.exponents[main_var] == *d { - let scale_factor = g.field.inv(t.coefficient); + let scale_factor = g.ring.inv(t.coefficient); g = g.mul_coeff(scale_factor); found = true; break; @@ -898,7 +898,7 @@ impl MultivariatePolynomial { } } - rhs.push(a.field.zero()); + rhs.push(a.ring.zero()); } } @@ -923,21 +923,21 @@ impl MultivariatePolynomial { let row_eval_first = &row_sample_values[shape_map[0]]; // assume first constant is 1, which will form the rhs of our equation - let actual_rhs = a.field.mul( + let actual_rhs = a.ring.mul( rhs_sec, - &a.field.pow(&row_eval_first[0], sample_index as u64 + 1), + &a.ring.pow(&row_eval_first[0], sample_index as u64 + 1), ); for aa in row_eval_sec { - gfm.push(a.field.pow(aa, sample_index as u64 + 1)); + gfm.push(a.ring.pow(aa, sample_index as u64 + 1)); } // place the scaling term variables at the end for aa in &row_eval_first[1..] { gfm.push( - a.field.neg( - &a.field - .mul(rhs_sec, &a.field.pow(aa, sample_index as u64 + 1)), + a.ring.neg( + &a.ring + .mul(rhs_sec, &a.ring.pow(aa, sample_index as u64 + 1)), ), ); } @@ -949,7 +949,7 @@ impl MultivariatePolynomial { // that yielded underdetermined systems for extra_relations in &scaling_var_relations { for _ in 0..vars_second { - gfm.push(a.field.zero()); + gfm.push(a.ring.zero()); } for v in &extra_relations[..vars_scale] { @@ -962,10 +962,10 @@ impl MultivariatePolynomial { gfm, rows as u32, samples_needed as u32, - a.field.clone(), + a.ring.clone(), ) .unwrap(); - let rhs = Matrix::new_vec(new_rhs, a.field.clone()); + let rhs = Matrix::new_vec(new_rhs, a.ring.clone()); match m.solve(&rhs) { Ok(r) => { @@ -1039,21 +1039,21 @@ impl MultivariatePolynomial { for sample_index in 0..max_terms { let row_eval_first = &row_sample_values[shape_map[0]]; let mut scaling_factor = - a.field.pow(&row_eval_first[0], sample_index as u64 + 1); // coeff eval is 1 + a.ring.pow(&row_eval_first[0], sample_index as u64 + 1); // coeff eval is 1 for (exp_eval, coeff_eval) in row_sample_values[shape_map[0]][1..].iter().zip(&r) { - a.field.add_mul_assign( + a.ring.add_mul_assign( &mut scaling_factor, coeff_eval, - &a.field.pow(exp_eval, sample_index as u64 + 1), + &a.ring.pow(exp_eval, sample_index as u64 + 1), ); } debug!( "Scaling fac {}: {}", sample_index, - a.field.printer(&scaling_factor) + a.ring.printer(&scaling_factor) ); lcoeff_cache.push(scaling_factor); } @@ -1061,7 +1061,7 @@ impl MultivariatePolynomial { for ((c, _), rhs) in shape.iter().zip(&mut samples) { rhs.truncate(c.nterms()); // drop unneeded samples for (r, scale) in rhs.iter_mut().zip(&lcoeff_cache) { - a.field.mul_assign(r, scale); + a.ring.mul_assign(r, scale); } } } else { @@ -1113,7 +1113,7 @@ impl, E: Exponent> MultivariatePolynomial { debug!("Content in last variable is not 1, but {}", c); // TODO: we assume that a content of -1 is also allowed // like in the special case gcd_(-x0*x1,-x0-x0*x1) - if c.nterms() != 1 || c.coefficients[0] != a.field.neg(&a.field.one()) { + if c.nterms() != 1 || c.coefficients[0] != a.ring.neg(&a.ring.one()) { return None; } } @@ -1138,26 +1138,26 @@ impl, E: Exponent> MultivariatePolynomial { } failure_count += 1; - if !a.field.size().is_zero() && failure_count * 2 > a.field.size() { + if !a.ring.size().is_zero() && failure_count * 2 > a.ring.size() { debug!("Cannot find unique sampling points: prime field is likely too small"); return None; } let mut sample_fail_count = 0i64; let v = loop { - let r = a.field.sample(&mut rng, (1, MAX_RNG_PREFACTOR as i64)); + let r = a.ring.sample(&mut rng, (1, MAX_RNG_PREFACTOR as i64)); if !gamma.replace(lastvar, &r).is_zero() { break r; } sample_fail_count += 1; - if !a.field.size().is_zero() && sample_fail_count * 2 > a.field.size() { + if !a.ring.size().is_zero() && sample_fail_count * 2 > a.ring.size() { debug!("Cannot find unique sampling points: prime field is likely too small"); continue 'newfirstnum; } }; - debug!("Chosen variable: {}", a.field.printer(&v)); + debug!("Chosen variable: {}", a.ring.printer(&v)); let av = a.replace(lastvar, &v); let bv = b.replace(lastvar, &v); @@ -1189,7 +1189,7 @@ impl, E: Exponent> MultivariatePolynomial { debug!( "GCD shape suggestion for sample point {} and gamma {}: {}", - a.field.printer(&v), + a.ring.printer(&v), gamma, gv ); @@ -1226,7 +1226,7 @@ impl, E: Exponent> MultivariatePolynomial { let mut gseq = vec![gv.clone().mul_coeff( gamma - .field + .ring .div(&gamma.replace(lastvar, &v).coefficients[0], &lc), )]; let mut vseq = vec![v]; @@ -1241,7 +1241,7 @@ impl, E: Exponent> MultivariatePolynomial { } let v = loop { - let v = a.field.sample(&mut rng, (1, MAX_RNG_PREFACTOR as i64)); + let v = a.ring.sample(&mut rng, (1, MAX_RNG_PREFACTOR as i64)); if !gamma.replace(lastvar, &v).is_zero() { // we need unique sampling points if !vseq.contains(&v) { @@ -1250,7 +1250,7 @@ impl, E: Exponent> MultivariatePolynomial { } sample_fail_count += 1; - if !a.field.size().is_zero() && sample_fail_count * 2 > a.field.size() { + if !a.ring.size().is_zero() && sample_fail_count * 2 > a.ring.size() { debug!( "Cannot find unique sampling points: prime field is likely too small" ); @@ -1258,7 +1258,7 @@ impl, E: Exponent> MultivariatePolynomial { } }; - debug!("Chosen sample: {}", a.field.printer(&v)); + debug!("Chosen sample: {}", a.ring.printer(&v)); let av = a.replace(lastvar, &v); let bv = b.replace(lastvar, &v); @@ -1303,7 +1303,7 @@ impl, E: Exponent> MultivariatePolynomial { debug!("Bad current image"); sample_fail_count += 1; - if !a.field.size().is_zero() && sample_fail_count * 2 > a.field.size() { + if !a.ring.size().is_zero() && sample_fail_count * 2 > a.ring.size() { debug!("Too many bad current images: prime field is likely too small"); continue 'newfirstnum; } @@ -1317,7 +1317,7 @@ impl, E: Exponent> MultivariatePolynomial { gseq.push( gv.clone().mul_coeff( gamma - .field + .ring .div(&gamma.replace(lastvar, &v).coefficients[0], &lc), ), ); @@ -1341,7 +1341,7 @@ impl, E: Exponent> MultivariatePolynomial { let mut cache = (0..a.nvars()) .map(|i| { vec![ - a.field.zero(); + a.ring.zero(); min( max(a.degree(i), b.degree(i)).to_u32() as usize + 1, POW_CACHE_SIZE @@ -1353,7 +1353,7 @@ impl, E: Exponent> MultivariatePolynomial { let r: Vec<_> = vars .iter() .skip(1) - .map(|i| (*i, a.field.sample(&mut rng, (1, MAX_RNG_PREFACTOR as i64)))) + .map(|i| (*i, a.ring.sample(&mut rng, (1, MAX_RNG_PREFACTOR as i64)))) .collect(); let g1 = gc.replace_all_except(vars[0], &r, &mut cache); @@ -1513,8 +1513,8 @@ impl, E: Exponent> MultivariatePolynomial< if self.is_constant() { let mut gcd = self.coefficients[0].clone(); for c in &b.coefficients { - gcd = self.field.gcd(&gcd, c); - if R::one_is_gcd_unit() && self.field.is_one(&gcd) { + gcd = self.ring.gcd(&gcd, c); + if R::one_is_gcd_unit() && self.ring.is_one(&gcd) { break; } } @@ -1524,8 +1524,8 @@ impl, E: Exponent> MultivariatePolynomial< if b.is_constant() { let mut gcd = b.coefficients[0].clone(); for c in &self.coefficients { - gcd = self.field.gcd(&gcd, c); - if R::one_is_gcd_unit() && self.field.is_one(&gcd) { + gcd = self.ring.gcd(&gcd, c); + if R::one_is_gcd_unit() && self.ring.is_one(&gcd) { break; } } @@ -1538,7 +1538,6 @@ impl, E: Exponent> MultivariatePolynomial< /// Compute the gcd of two multivariate polynomials. #[instrument(skip_all)] pub fn gcd(&self, b: &MultivariatePolynomial) -> MultivariatePolynomial { - debug_assert_eq!(self.nvars(), b.nvars()); debug!("gcd of {} and {}", self, b); if let Some(g) = self.simple_gcd(b) { @@ -1550,6 +1549,10 @@ impl, E: Exponent> MultivariatePolynomial< let mut a = Cow::Borrowed(self); let mut b = Cow::Borrowed(b); + if self.variables != b.variables { + a.to_mut().unify_variables(b.to_mut()); + } + // determine the maximum shared power of every variable let mut shared_degree: SmallVec<[E; INLINED_EXPONENTS]> = a.exponents(0).into(); for p in [&a, &b] { @@ -1580,7 +1583,7 @@ impl, E: Exponent> MultivariatePolynomial< let mut base_degree: SmallVec<[Option; INLINED_EXPONENTS]> = smallvec![None; a.nvars()]; if let Some(g) = MultivariatePolynomial::simple_gcd(&a, &b) { - return rescale_gcd(g, &shared_degree, &base_degree, &a.constant(a.field.one())); + return rescale_gcd(g, &shared_degree, &base_degree, &a.constant(a.ring.one())); } // check if the polynomial are functions of x^n, n > 1 @@ -1661,7 +1664,7 @@ impl, E: Exponent> MultivariatePolynomial< gcd.0, &shared_degree, &base_degree, - &a.constant(a.field.one()), + &a.constant(a.ring.one()), ); } @@ -1738,7 +1741,7 @@ impl, E: Exponent> MultivariatePolynomial< cont, &shared_degree, &base_degree, - &p1.constant(p1.field.one()), + &p1.constant(p1.ring.one()), ); } } @@ -1790,13 +1793,13 @@ impl, E: Exponent> MultivariatePolynomial< // get the integer content for univariate polynomials let uca = a.content(); let ucb = b.content(); - let content = a.field.gcd(&a.content(), &b.content()); + let content = a.ring.gcd(&a.content(), &b.content()); let p = a.zero_with_capacity(1); - if !a.field.is_one(&uca) { + if !a.ring.is_one(&uca) { a = Cow::Owned(a.into_owned().div_coeff(&uca)); } - if !a.field.is_one(&ucb) { + if !a.ring.is_one(&ucb) { b = Cow::Owned(b.into_owned().div_coeff(&ucb)); } @@ -1903,14 +1906,14 @@ impl MultivariatePolynomial { debug!("a={}; b={}", self, b); // do integer GCD - let content_gcd = self.field.gcd(&self.content(), &b.content()); + let content_gcd = self.ring.gcd(&self.content(), &b.content()); debug!("content={}", content_gcd); let mut a = Cow::Borrowed(self); let mut b = Cow::Borrowed(b); - if !a.field.is_one(&content_gcd) { + if !a.ring.is_one(&content_gcd) { a = Cow::Owned(a.into_owned().div_coeff(&content_gcd)); b = Cow::Owned(b.into_owned().div_coeff(&content_gcd)); } @@ -2040,11 +2043,11 @@ impl MultivariatePolynomial { if let Some(n) = f.iter().find(|x| x.is_constant()) { let mut gcd = n.content(); for x in f.iter() { - if x.field.is_one(&gcd) { + if x.ring.is_one(&gcd) { break; } - gcd = x.field.gcd(&gcd, &x.content()); + gcd = x.ring.gcd(&gcd, &x.content()); } return n.constant(gcd); } @@ -2092,7 +2095,7 @@ impl MultivariatePolynomial { f.retain(|x| { if x.divides(&gcd).is_some() { - content_gcd = gcd.field.gcd(&content_gcd, &x.content()); + content_gcd = gcd.ring.gcd(&content_gcd, &x.content()); false } else { true @@ -2145,7 +2148,7 @@ impl MultivariatePolynomial { // compute scaling factor in Z let gamma = self - .field + .ring .gcd(&self.lcoeff_varorder(vars), &b.lcoeff_varorder(vars)); debug!("gamma {}", gamma); @@ -2226,7 +2229,7 @@ impl MultivariatePolynomial { } let gpc = gp.lcoeff_varorder(vars); - let lcoeff_factor = gp.field.div(&gammap, &gpc); + let lcoeff_factor = gp.ring.div(&gammap, &gpc); // construct the gcd suggestion in Z let mut gm = self.zero_with_capacity(gp.nterms()); @@ -2235,8 +2238,8 @@ impl MultivariatePolynomial { .coefficients .iter() .map(|x| { - gp.field - .to_symmetric_integer(&gp.field.mul(x, &lcoeff_factor)) + gp.ring + .to_symmetric_integer(&gp.ring.mul(x, &lcoeff_factor)) }) .collect(); @@ -2351,8 +2354,8 @@ impl MultivariatePolynomial { // scale the new image let gpc = gp.lcoeff_varorder(vars); - gp = gp.mul_coeff(ap.field.div(&gammap, &gpc)); - debug!("gp: {} mod {}", gp, gp.field.get_prime()); + gp = gp.mul_coeff(ap.ring.div(&gammap, &gpc)); + debug!("gp: {} mod {}", gp, gp.ring.get_prime()); // use chinese remainder theorem to merge coefficients and map back to Z // terms could be missing in gp, but not in gm (TODO: check this?) @@ -2362,26 +2365,26 @@ impl MultivariatePolynomial { gpi += 1; gp.coefficients[gpi - 1] } else { - ap.field.zero() + ap.ring.zero() }; let gmc = &mut gm.coefficients[t]; let coeff = if gmc.is_negative() { - self.field.add(gmc, &m) + self.ring.add(gmc, &m) } else { gmc.clone() }; *gmc = Integer::chinese_remainder( coeff, - Integer::from_finite_field(&gp.field, gpc), + Integer::from_finite_field(&gp.ring, gpc), m.clone(), - Integer::from_prime(&gp.field), + Integer::from_prime(&gp.ring), ); } - self.field - .mul_assign(&mut m, &Integer::from_prime(&gp.field)); + self.ring + .mul_assign(&mut m, &Integer::from_prime(&gp.ring)); debug!("gm: {} from ring {}", gm, m); } @@ -2554,10 +2557,10 @@ impl PolynomialGCD for RationalField { tight_bounds: &mut [E], ) -> MultivariatePolynomial { // remove the content so that the polynomials have integer coefficients - let content = a.field.gcd(&a.content(), &b.content()); + let content = a.ring.gcd(&a.content(), &b.content()); - let a_int = a.map_coeff(|c| a.field.div(c, &content).numerator(), Z); - let b_int = b.map_coeff(|c| b.field.div(c, &content).numerator(), Z); + let a_int = a.map_coeff(|c| a.ring.div(c, &content).numerator(), Z); + let b_int = b.map_coeff(|c| b.ring.div(c, &content).numerator(), Z); MultivariatePolynomial::gcd_zippel::(&a_int, &b_int, vars, bounds, tight_bounds) .map_coeff(|c| c.to_rational(), Q) @@ -2635,11 +2638,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 - 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 field = a.ring.upgrade(a.ring.get_extension_degree() as usize + 1); + let ag = a.map_coeff(|c| a.ring.upgrade_element(c, &field), field.clone()); + let bg = b.map_coeff(|c| a.ring.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()) + g.map_coeff(|c| a.ring.downgrade_element(c), a.ring.clone()) } } } @@ -2691,9 +2694,9 @@ impl PolynomialGCD for AlgebraicExtension { bounds: &mut [E], tight_bounds: &mut [E], ) -> MultivariatePolynomial { - let content = a.field.poly().content().inv(); + let content = a.ring.poly().content().inv(); let a_integer = - AlgebraicExtension::new(a.field.poly().map_coeff(|c| (c * &content).numerator(), Z)); + AlgebraicExtension::new(a.ring.poly().map_coeff(|c| (c * &content).numerator(), Z)); let a_lcoeff = a_integer.poly().lcoeff(); debug!("Zippel gcd of {} and {} % {}", a, b, a_integer); @@ -2721,7 +2724,7 @@ impl PolynomialGCD for AlgebraicExtension { let mut p = &primes[pi]; let mut finite_field = Zp::new(*p); - let mut algebraic_field_ff = a.field.to_finite_field(&finite_field); + let mut algebraic_field_ff = a.ring.to_finite_field(&finite_field); let a_lcoeff_p = a_lcoeff.to_finite_field(&finite_field); @@ -2789,7 +2792,7 @@ impl PolynomialGCD for AlgebraicExtension { } let gpc = gp.lcoeff_varorder(vars); - let lcoeff_factor = gp.field.inv(&gpc); + let lcoeff_factor = gp.ring.inv(&gpc); // construct the gcd suggestion in Z // contrary to the integer case, we do not know the leading coefficient in Z @@ -2803,7 +2806,7 @@ impl PolynomialGCD for AlgebraicExtension { .iter() .map(|x| { a_integer.to_element( - gp.field + gp.ring .mul(x, &lcoeff_factor) .poly .map_coeff(|c| finite_field.to_symmetric_integer(c), Z), @@ -2831,7 +2834,7 @@ impl PolynomialGCD for AlgebraicExtension { p = &primes[pi]; finite_field = Zp::new(*p); - algebraic_field_ff = a.field.to_finite_field(&finite_field); + algebraic_field_ff = a.ring.to_finite_field(&finite_field); let a_lcoeff_p = a_lcoeff.to_finite_field(&finite_field); @@ -2911,8 +2914,8 @@ impl PolynomialGCD for AlgebraicExtension { // scale the new image let gpc = gp.lcoeff_varorder(vars); - gp = gp.mul_coeff(ap.field.inv(&gpc)); - debug!("gp: {} mod {}", gp, gp.field); + gp = gp.mul_coeff(ap.ring.inv(&gpc)); + debug!("gp: {} mod {}", gp, gp.ring); // use chinese remainder theorem to merge coefficients and map back to Z // terms could be missing in gp, but not in gm (TODO: check this?) @@ -2922,7 +2925,7 @@ impl PolynomialGCD for AlgebraicExtension { gpi += 1; gp.coefficients[gpi - 1].clone() } else { - ap.field.zero() + ap.ring.zero() }; let gmc_a = &mut gm.coefficients[t]; @@ -2930,7 +2933,7 @@ impl PolynomialGCD for AlgebraicExtension { // apply CRT to each integer coefficient in the algebraic number ring let mut gpc_pos = 0; let mut gmc_pos = 0; - for i in 0..a.field.poly().degree(0) { + for i in 0..a.ring.poly().degree(0) { let gpc = if gpc_pos < gpc.poly.nterms() && i == gpc.poly.exponents(gpc_pos)[0] { gpc_pos += 1; @@ -2986,7 +2989,7 @@ impl PolynomialGCD for AlgebraicExtension { let mut gc = a.zero(); for c in &gm.coefficients { - let mut nc = a.field.poly().zero(); + let mut nc = a.ring.poly().zero(); for aa in &c.poly.coefficients { match Rational::maximal_quotient_reconstruction(aa, &m, None) { @@ -3000,7 +3003,7 @@ impl PolynomialGCD for AlgebraicExtension { } nc.exponents.clone_from(&c.poly.exponents); - gc.coefficients.push(a.field.to_element(nc)); + gc.coefficients.push(a.ring.to_element(nc)); } gc.exponents.clone_from(&gm.exponents); @@ -3026,7 +3029,7 @@ impl PolynomialGCD for AlgebraicExtension { let mut i = 0; loop { let f = Zp::new(LARGE_U32_PRIMES[i]); - let algebraic_field_ff = a.field.to_finite_field(&f); + let algebraic_field_ff = a.ring.to_finite_field(&f); let ap = a.map_coeff(|c| c.to_finite_field(&f), algebraic_field_ff.clone()); let bp = b.map_coeff(|c| c.to_finite_field(&f), algebraic_field_ff.clone()); if ap.nterms() > 0 diff --git a/src/poly/groebner.rs b/src/poly/groebner.rs index 2533a999..95c242b1 100644 --- a/src/poly/groebner.rs +++ b/src/poly/groebner.rs @@ -121,7 +121,7 @@ impl GroebnerBasis GroebnerBasis { *e = *e1 - *e2; } - let ratio = g.field.div(r.max_coeff(), g.max_coeff()); + let ratio = g.ring.div(r.max_coeff(), g.max_coeff()); r = r - g.clone().mul_exp(&monom).mul_coeff(ratio); if r.is_zero() { @@ -507,7 +507,7 @@ impl GroebnerBasis { lead_reduced.swap(0, i); let h = Self::reduce(&lead_reduced[0], &lead_reduced[1..]); if !h.is_zero() { - let i = h.field.inv(h.max_coeff()); + let i = h.ring.inv(h.max_coeff()); basis.push(h.mul_coeff(i)); } } @@ -545,11 +545,11 @@ impl GroebnerBasis { let new_f1 = p1 .clone() .mul_exp(&extra_factor_f1) - .mul_coeff(p1.field.div(p2.max_coeff(), p1.max_coeff())); + .mul_coeff(p1.ring.div(p2.max_coeff(), p1.max_coeff())); let new_f2 = p2 .clone() .mul_exp(&extra_factor_f2) - .mul_coeff(p1.field.div(p1.max_coeff(), p2.max_coeff())); + .mul_coeff(p1.ring.div(p1.max_coeff(), p2.max_coeff())); let s = new_f1 - new_f2; diff --git a/src/poly/polynomial.rs b/src/poly/polynomial.rs index c2e4bcd2..06255e61 100755 --- a/src/poly/polynomial.rs +++ b/src/poly/polynomial.rs @@ -25,23 +25,20 @@ thread_local! { static DENSE_MUL_BUFFER: Cell> = const { Cell::new(Vec: #[derive(Clone, PartialEq, Eq, Hash, Debug)] pub struct PolynomialRing { pub(crate) ring: R, - pub(crate) variables: Arc>, _phantom_exp: PhantomData, } impl PolynomialRing { - pub fn new(coeff_ring: R, var_map: Arc>) -> PolynomialRing { + pub fn new(coeff_ring: R) -> PolynomialRing { PolynomialRing { ring: coeff_ring, - variables: var_map, _phantom_exp: PhantomData, } } - pub fn new_from_poly(poly: &MultivariatePolynomial) -> PolynomialRing { + pub fn from_poly(poly: &MultivariatePolynomial) -> PolynomialRing { PolynomialRing { - ring: poly.field.clone(), - variables: poly.variables.clone(), + ring: poly.ring.clone(), _phantom_exp: PhantomData, } } @@ -108,7 +105,7 @@ impl Ring for PolynomialRing { #[inline] fn zero(&self) -> Self::Element { - MultivariatePolynomial::new(&self.ring, None, self.variables.clone()) + MultivariatePolynomial::new(&self.ring, None, Arc::new(vec![])) } #[inline] @@ -198,7 +195,7 @@ impl, E: Exponent> EuclideanDomain for Pol } } -/// Multivariate polynomial with a sparse degree and variable dense representation. +/// Multivariate polynomial with a sparse degree and dense variable representation. #[derive(Clone)] pub struct MultivariatePolynomial { // Data format: the i-th monomial is stored as coefficients[i] and @@ -206,7 +203,8 @@ pub struct MultivariatePolynomial, pub exponents: Vec, - pub field: F, + /// The coefficient ring. + pub ring: F, pub variables: Arc>, pub(crate) _phantom: PhantomData, } @@ -216,23 +214,51 @@ impl MultivariatePolynomial { /// prefer to create new polynomials from existing ones, so that the /// variable map and field are inherited. #[inline] - pub fn new(field: &F, cap: Option, variables: Arc>) -> Self { + pub fn new(ring: &F, cap: Option, variables: Arc>) -> Self { Self { coefficients: Vec::with_capacity(cap.unwrap_or(0)), exponents: Vec::with_capacity(cap.unwrap_or(0) * variables.len()), - field: field.clone(), + ring: ring.clone(), variables, _phantom: PhantomData, } } + /// Constructs a zero polynomial. Instead of using this constructor, + /// prefer to create new polynomials from existing ones, so that the + /// variable map is inherited. + #[inline] + pub fn new_zero(ring: &F) -> Self { + Self { + coefficients: vec![], + exponents: vec![], + ring: ring.clone(), + variables: Arc::new(vec![]), + _phantom: PhantomData, + } + } + + /// Constructs a polynomial that is one. Instead of using this constructor, + /// prefer to create new polynomials from existing ones, so that the + /// variable map is inherited. + #[inline] + pub fn new_one(ring: &F) -> Self { + Self { + coefficients: vec![ring.one()], + exponents: vec![], + ring: ring.clone(), + variables: Arc::new(vec![]), + _phantom: PhantomData, + } + } + /// Constructs a zero polynomial, inheriting the field and variable map from `self`. #[inline] pub fn zero(&self) -> Self { Self { coefficients: vec![], exponents: vec![], - field: self.field.clone(), + ring: self.ring.clone(), variables: self.variables.clone(), _phantom: PhantomData, } @@ -245,7 +271,7 @@ impl MultivariatePolynomial { Self { coefficients: Vec::with_capacity(cap), exponents: Vec::with_capacity(cap * self.nvars()), - field: self.field.clone(), + ring: self.ring.clone(), variables: self.variables.clone(), _phantom: PhantomData, } @@ -262,7 +288,7 @@ impl MultivariatePolynomial { Self { coefficients: vec![coeff], exponents: vec![E::zero(); self.nvars()], - field: self.field.clone(), + ring: self.ring.clone(), variables: self.variables.clone(), _phantom: PhantomData, } @@ -272,9 +298,9 @@ impl MultivariatePolynomial { #[inline] pub fn one(&self) -> Self { Self { - coefficients: vec![self.field.one()], + coefficients: vec![self.ring.one()], exponents: vec![E::zero(); self.nvars()], - field: self.field.clone(), + ring: self.ring.clone(), variables: self.variables.clone(), _phantom: PhantomData, } @@ -292,7 +318,7 @@ impl MultivariatePolynomial { Self { coefficients: vec![coeff], exponents, - field: self.field.clone(), + ring: self.ring.clone(), variables: self.variables.clone(), _phantom: PhantomData, } @@ -304,7 +330,7 @@ impl MultivariatePolynomial { if let Some(pos) = self.variables.iter().position(|v| v == var) { let mut exp = vec![E::zero(); self.nvars()]; exp[pos] = E::one(); - Ok(self.monomial(self.field.one(), exp)) + Ok(self.monomial(self.ring.one(), exp)) } else { Err("Variable not found") } @@ -335,7 +361,7 @@ impl MultivariatePolynomial { #[inline] pub fn is_one(&self) -> bool { self.nterms() == 1 - && self.field.is_one(&self.coefficients[0]) + && self.ring.is_one(&self.coefficients[0]) && self.exponents.iter().all(|x| x.is_zero()) } @@ -368,7 +394,7 @@ impl MultivariatePolynomial { #[inline] pub fn get_constant(&self) -> F::Element { if self.is_zero() || !self.exponents(0).iter().all(|e| e.is_zero()) { - return self.field.zero(); + return self.ring.zero(); } self.coefficients[0].clone() @@ -489,7 +515,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.ring, 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 { @@ -598,7 +624,7 @@ impl MultivariatePolynomial { let nterms = self.nterms(); if nterms > 0 && exponents == self.last_exponents() { - self.field + self.ring .add_assign(&mut self.coefficients[nterms - 1], &coefficient); if F::is_zero(&self.coefficients[nterms - 1]) { @@ -648,7 +674,7 @@ impl MultivariatePolynomial { match c { Ordering::Equal => { // Add the two coefficients. - self.field + self.ring .add_assign(&mut self.coefficients[m], &coefficient); if F::is_zero(&self.coefficients[m]) { // The coefficient becomes zero. Remove this monomial. @@ -696,7 +722,7 @@ impl MultivariatePolynomial { exp.copy_from_slice(x.exponents); let pow = exp[var].to_u32() as u64; exp[var] = exp[var] - E::one(); - res.append_monomial(self.field.mul(x.coefficient, &self.field.nth(pow)), &exp); + res.append_monomial(self.ring.mul(x.coefficient, &self.ring.nth(pow)), &exp); } } @@ -738,7 +764,7 @@ impl Display for MultivariateP impl PartialEq for MultivariatePolynomial { #[inline] fn eq(&self, other: &Self) -> bool { - if self.nvars() != other.nvars() { + if self.variables != other.variables { if self.is_constant() != other.is_constant() { return false; } @@ -755,10 +781,7 @@ impl PartialEq for MultivariatePolynomia return self.coefficients[0] == other.coefficients[0]; } - // TODO: what is expected here? - unimplemented!( - "Cannot compare non-constant polynomials with different variable maps yet" - ); + return false; } if self.nterms() != other.nterms() { return false; @@ -771,7 +794,10 @@ impl std::hash::Hash for MultivariatePol fn hash(&self, state: &mut H) { self.coefficients.hash(state); self.exponents.hash(state); - self.variables.hash(state); + + if !self.is_constant() { + self.variables.hash(state); + } } } @@ -780,6 +806,7 @@ impl Eq for MultivariatePolynomial InternalOrdering for MultivariatePolynomial { /// An ordering of polynomials that has no intuitive meaning. fn internal_cmp(&self, other: &Self) -> Ordering { + // TODO: what about different variables? Ord::cmp(&self.exponents, &other.exponents) .then_with(|| self.coefficients.internal_cmp(&other.coefficients)) } @@ -789,7 +816,7 @@ impl Add for MultivariatePolynomial Self::Output { - debug_assert_eq!(self.field, other.field); + assert_eq!(self.ring, other.ring); self.unify_variables(&mut other); @@ -802,7 +829,7 @@ impl Add for MultivariatePolynomial = vec![E::zero(); self.nvars() * (self.nterms() + other.nterms())]; let mut new_nterms = 0; @@ -834,7 +861,7 @@ impl Add for MultivariatePolynomial { - self.field + self.ring .add_assign(&mut self.coefficients[i], &other.coefficients[j]); if !F::is_zero(&self.coefficients[i]) { insert_monomial!(self, i); @@ -861,7 +888,7 @@ impl Add for MultivariatePolynomial Add<&'a MultivariatePolynom type Output = MultivariatePolynomial; fn add(self, other: &MultivariatePolynomial) -> Self::Output { - debug_assert_eq!(self.field, other.field); - - if self.variables != other.variables { - let mut c1 = self.clone(); - let mut c2 = other.clone(); - c1.unify_variables(&mut c2); - return c1 + c2; - } + assert_eq!(self.ring, other.ring); if self.is_zero() { return other.clone(); @@ -890,8 +910,15 @@ impl<'a, 'b, F: Ring, E: Exponent, O: MonomialOrder> Add<&'a MultivariatePolynom return self.clone(); } + if self.variables != other.variables { + let mut c1 = self.clone(); + let mut c2 = other.clone(); + c1.unify_variables(&mut c2); + return c1 + c2; + } + // Merge the two polynomials, which are assumed to be already sorted. - let mut new_coefficients = vec![self.field.zero(); self.nterms() + other.nterms()]; + let mut new_coefficients = vec![self.ring.zero(); self.nterms() + other.nterms()]; let mut new_exponents: Vec = vec![E::zero(); self.nvars() * (self.nterms() + other.nterms())]; let mut new_nterms = 0; @@ -919,9 +946,7 @@ impl<'a, 'b, F: Ring, E: Exponent, O: MonomialOrder> Add<&'a MultivariatePolynom j += 1; } Ordering::Equal => { - let coeff = self - .field - .add(&self.coefficients[i], &other.coefficients[j]); + let coeff = self.ring.add(&self.coefficients[i], &other.coefficients[j]); if !F::is_zero(&coeff) { new_coefficients[new_nterms] = coeff; new_exponents[new_nterms * self.nvars()..(new_nterms + 1) * self.nvars()] @@ -950,7 +975,7 @@ impl<'a, 'b, F: Ring, E: Exponent, O: MonomialOrder> Add<&'a MultivariatePolynom MultivariatePolynomial { coefficients: new_coefficients, exponents: new_exponents, - field: self.field.clone(), + ring: self.ring.clone(), variables: self.variables.clone(), _phantom: PhantomData, } @@ -980,7 +1005,7 @@ impl Neg for MultivariatePolynomial Self::Output { // Negate coefficients of all terms. for c in &mut self.coefficients { - *c = self.field.neg(c); + *c = self.ring.neg(c); } self } @@ -993,10 +1018,26 @@ impl<'a, 'b, F: Ring, E: Exponent> Mul<&'a MultivariatePolynomial) -> Self::Output { + assert_eq!(self.ring, rhs.ring); + if self.nterms() == 0 || rhs.nterms() == 0 { return self.zero(); } + if self.is_constant() { + return rhs.clone().mul_coeff(self.coefficients[0].clone()); + } + if rhs.is_constant() { + return self.clone().mul_coeff(rhs.coefficients[0].clone()); + } + + if self.variables != rhs.variables { + let mut c1 = self.clone(); + let mut c2 = rhs.clone(); + c1.unify_variables(&mut c2); + return c1.mul(&c2); + } + if self.nterms() == 1 { return rhs .clone() @@ -1072,7 +1113,7 @@ impl MultivariatePolynomial { MultivariatePolynomial { coefficients, exponents, - field: self.field.clone(), + ring: self.ring.clone(), variables: self.variables.clone(), _phantom: PhantomData, } @@ -1081,7 +1122,7 @@ impl MultivariatePolynomial { /// Multiply every coefficient with `other`. pub fn mul_coeff(mut self, other: F::Element) -> Self { for c in &mut self.coefficients { - self.field.mul_assign(c, &other); + self.ring.mul_assign(c, &other); } for i in (0..self.nterms()).rev() { @@ -1115,7 +1156,7 @@ impl MultivariatePolynomial { MultivariatePolynomial { coefficients, exponents, - field, + ring: field, variables: self.variables.clone(), _phantom: PhantomData, } @@ -1199,7 +1240,7 @@ impl MultivariatePolynomial { /// Get the leading coefficient. pub fn lcoeff(&self) -> F::Element { if self.is_zero() { - return self.field.zero(); + return self.ring.zero(); } self.coefficients.last().unwrap().clone() } @@ -1225,7 +1266,7 @@ impl MultivariatePolynomial { } let mut highest = vec![E::zero(); self.nvars()]; - let mut highestc = &self.field.zero(); + let mut highestc = &self.ring.zero(); 'nextmon: for m in self.into_iter() { let mut more = false; @@ -1431,7 +1472,7 @@ impl MultivariatePolynomial { 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()); + MultivariatePolynomial::new(&self.ring, self.nterms().into(), self.variables.clone()); for i in indices { res.append_monomial( @@ -1460,9 +1501,9 @@ impl MultivariatePolynomial { continue; } - let c = self.field.mul( + let c = self.ring.mul( t.coefficient, - &self.field.pow(v, t.exponents[n].to_u32() as u64), + &self.ring.pow(v, t.exponents[n].to_u32() as u64), ); e.copy_from_slice(t.exponents); @@ -1486,9 +1527,9 @@ impl MultivariatePolynomial { continue; } - let c = self.field.mul( + let c = self.ring.mul( t.coefficient, - &self.field.pow(v, t.exponents[n].to_u32() as u64), + &self.ring.pow(v, t.exponents[n].to_u32() as u64), ); if F::is_zero(&c) { @@ -1503,7 +1544,7 @@ impl MultivariatePolynomial { res.exponents.extend_from_slice(&e); } else { let l = res.coefficients.last_mut().unwrap(); - self.field.add_assign(l, &c); + self.ring.add_assign(l, &c); if F::is_zero(l) { res.coefficients.pop(); @@ -1518,7 +1559,7 @@ impl MultivariatePolynomial { /// Replace a variable `n` in the polynomial by an element from /// the ring `v`. pub fn replace_all(&self, r: &[F::Element]) -> F::Element { - let mut res = self.field.zero(); + let mut res = self.ring.zero(); // TODO: cache power taking? for t in self { @@ -1526,12 +1567,12 @@ impl MultivariatePolynomial { for (i, v) in r.iter().zip(t.exponents) { if v != &E::zero() { - self.field - .mul_assign(&mut c, &self.field.pow(i, v.to_u32() as u64)); + self.ring + .mul_assign(&mut c, &self.ring.pow(i, v.to_u32() as u64)); } } - self.field.add_assign(&mut res, &c); + self.ring.add_assign(&mut res, &c); } res @@ -1581,18 +1622,18 @@ impl MultivariatePolynomial { if p > 0 { if p < cache[*n].len() { if F::is_zero(&cache[*n][p]) { - cache[*n][p] = self.field.pow(vv, p as u64); + cache[*n][p] = self.ring.pow(vv, p as u64); } - self.field.mul_assign(&mut c, &cache[*n][p]); + self.ring.mul_assign(&mut c, &cache[*n][p]); } else { - self.field.mul_assign(&mut c, &self.field.pow(vv, p as u64)); + self.ring.mul_assign(&mut c, &self.ring.pow(vv, p as u64)); } } } tm.entry(t.exponents[v]) - .and_modify(|e| self.field.add_assign(e, &c)) + .and_modify(|e| self.ring.add_assign(e, &c)) .or_insert(c); } @@ -1614,7 +1655,7 @@ impl MultivariatePolynomial { } if self.is_constant() { - return self.constant(self.field.pow(&self.lcoeff(), pow as u64)); + return self.constant(self.ring.pow(&self.lcoeff(), pow as u64)); } let mut x = self.clone(); @@ -1636,7 +1677,7 @@ impl MultivariatePolynomial { let c = self.to_univariate_polynomial_list(var); let mut p = UnivariatePolynomial::new( - &PolynomialRing::new_from_poly(self), + &PolynomialRing::from_poly(self), None, Arc::new(self.variables[var].clone()), ); @@ -1645,7 +1686,7 @@ impl MultivariatePolynomial { return p; } - p.coefficients = vec![p.field.zero(); c.last().unwrap().1.to_u32() as usize + 1]; + p.coefficients = vec![self.zero(); c.last().unwrap().1.to_u32() as usize + 1]; for (q, e) in c { p.coefficients[e.to_u32() as usize] = q; } @@ -1655,7 +1696,7 @@ impl MultivariatePolynomial { pub fn to_univariate_from_univariate(&self, var: usize) -> UnivariatePolynomial { let mut p = - UnivariatePolynomial::new(&self.field, None, Arc::new(self.variables[var].clone())); + UnivariatePolynomial::new(&self.ring, None, Arc::new(self.variables[var].clone())); if self.is_zero() { return p; @@ -1787,7 +1828,7 @@ impl MultivariatePolynomial { MultivariatePolynomial::new(field, self.nterms().into(), self.variables.clone()); for (e, c) in polys { let mut c2 = MultivariatePolynomial::new( - &self.field, + &self.ring, c.nterms().into(), Arc::new(vec![self.variables.as_ref()[var_index].clone()]), ); @@ -1843,7 +1884,7 @@ impl MultivariatePolynomial { max = max.min(m); } - let mut coeffs = vec![self.field.zero(); max + 1]; + let mut coeffs = vec![self.ring.zero(); max + 1]; for x in self { for y in rhs { @@ -1852,7 +1893,7 @@ impl MultivariatePolynomial { continue; } - self.field + self.ring .add_mul_assign(&mut coeffs[pos as usize], x.coefficient, y.coefficient); } } @@ -1936,12 +1977,12 @@ impl MultivariatePolynomial { // check if we need to use a dense indexing array to save memory if total < 1000 { - let mut coeffs = vec![self.field.zero(); total]; + let mut coeffs = vec![self.ring.zero(); total]; for (c1, e1) in self.coefficients.iter().zip(&uni_exp_self) { for (c2, e2) in rhs.coefficients.iter().zip(&uni_exp_rhs) { let pos = e1.to_u32() as usize + e2.to_u32() as usize; - self.field.add_mul_assign(&mut coeffs[pos], c1, c2); + self.ring.add_mul_assign(&mut coeffs[pos], c1, c2); } } @@ -1966,10 +2007,10 @@ impl MultivariatePolynomial { for (c2, e2) in rhs.coefficients.iter().zip(&uni_exp_rhs) { let pos = e1.to_u32() as usize + e2.to_u32() as usize; if coeff_index[pos] == 0 { - coeffs.push(self.field.mul(c1, c2)); + coeffs.push(self.ring.mul(c1, c2)); coeff_index[pos] = coeffs.len() as u32; } else { - self.field.add_mul_assign( + self.ring.add_mul_assign( &mut coeffs[coeff_index[pos] as usize - 1], c1, c2, @@ -1982,7 +2023,7 @@ impl MultivariatePolynomial { if *c != 0 { from_uni_var(p as u32, &max_degs_rev, &mut exp); r.append_monomial( - std::mem::replace(&mut coeffs[*c as usize - 1], self.field.zero()), + std::mem::replace(&mut coeffs[*c as usize - 1], self.ring.zero()), &exp, ); *c = 0; @@ -2062,12 +2103,12 @@ impl MultivariatePolynomial { while !h.is_empty() { let cur_mon = h.pop().unwrap(); - let mut coefficient = self.field.zero(); + let mut coefficient = self.ring.zero(); let mut q = cache.remove(&cur_mon.0).unwrap(); for (i, j) in q.drain(..) { - self.field.add_mul_assign( + self.ring.add_mul_assign( &mut coefficient, &self.coefficients[i], &rhs.coefficients[j], @@ -2173,12 +2214,12 @@ impl MultivariatePolynomial { in_heap[0] = true; while let Some(cur_mon) = h.pop() { - let mut coefficient = self.field.zero(); + let mut coefficient = self.ring.zero(); let mut q = cache.remove(&cur_mon.0).unwrap(); for (i, j) in q.drain(..) { - self.field.add_mul_assign( + self.ring.add_mul_assign( &mut coefficient, &self.coefficients[i], &other.coefficients[j], @@ -2249,7 +2290,7 @@ impl MultivariatePolynomial { MultivariatePolynomial, MultivariatePolynomial, ) { - debug_assert_eq!(div.lcoeff(), self.field.one()); + debug_assert_eq!(div.lcoeff(), self.ring.one()); if self.is_zero() { return (self.clone(), self.clone()); } @@ -2278,7 +2319,7 @@ impl MultivariatePolynomial { break self.coefficients[dividendpos].clone(); } if dividendpos == 0 || self.exponents(dividendpos)[var] < pow { - break self.field.zero(); + break self.ring.zero(); } dividendpos -= 1; }; @@ -2293,7 +2334,7 @@ impl MultivariatePolynomial { } if div.exponents(bindex)[var] + q.exponents(qindex)[var] == pow { - self.field.sub_mul_assign( + self.ring.sub_mul_assign( &mut coeff, &div.coefficients[bindex], &q.coefficients[qindex], @@ -2383,16 +2424,16 @@ impl MultivariatePolynomial { /// Get the content from the coefficients. pub fn content(&self) -> F::Element { if self.coefficients.is_empty() { - return self.field.zero(); + return self.ring.zero(); } let mut c = self.coefficients.first().unwrap().clone(); for cc in self.coefficients.iter().skip(1) { // early return if possible (not possible for rationals) - if F::one_is_gcd_unit() && self.field.is_one(&c) { + if F::one_is_gcd_unit() && self.ring.is_one(&c) { break; } - c = self.field.gcd(&c, cc); + c = self.ring.gcd(&c, cc); } c } @@ -2400,7 +2441,7 @@ impl MultivariatePolynomial { /// Divide every coefficient with `other`. pub fn div_coeff(mut self, other: &F::Element) -> Self { for c in &mut self.coefficients { - let (quot, rem) = self.field.quot_rem(c, other); + let (quot, rem) = self.ring.quot_rem(c, other); debug_assert!(F::is_zero(&rem)); *c = quot; } @@ -2421,12 +2462,19 @@ impl MultivariatePolynomial { panic!("Cannot divide by 0 polynomial"); } + if self.variables != div.variables { + let mut c1 = self.clone(); + let mut c2 = div.clone(); + c1.unify_variables(&mut c2); + return self.divides(&c2); + } + if self.is_zero() { return Some(self.clone()); } // check if the leading coefficients divide - if !F::is_zero(&self.field.rem(&self.lcoeff(), &div.lcoeff())) { + if !F::is_zero(&self.ring.rem(&self.lcoeff(), &div.lcoeff())) { return None; } @@ -2434,29 +2482,29 @@ impl MultivariatePolynomial { return None; } - if self.field.characteristic().is_zero() { + if self.ring.characteristic().is_zero() { // test division of constant term (evaluation at x_i = 0) let c = div.get_constant(); if !F::is_zero(&c) - && !self.field.is_one(&c) - && !F::is_zero(&self.field.rem(&self.get_constant(), &c)) + && !self.ring.is_one(&c) + && !F::is_zero(&self.ring.rem(&self.get_constant(), &c)) { return None; } // test division at x_i = 1 - let mut num = self.field.zero(); + let mut num = self.ring.zero(); for c in &self.coefficients { - self.field.add_assign(&mut num, c); + self.ring.add_assign(&mut num, c); } - let mut den = self.field.zero(); + let mut den = self.ring.zero(); for c in &div.coefficients { - self.field.add_assign(&mut den, c); + self.ring.add_assign(&mut den, c); } if !F::is_zero(&den) - && !self.field.is_one(&den) - && !F::is_zero(&self.field.rem(&num, &den)) + && !self.ring.is_one(&den) + && !F::is_zero(&self.ring.rem(&num, &den)) { return None; } @@ -2496,13 +2544,20 @@ impl MultivariatePolynomial { return (self.clone(), self.zero()); } + if self.variables != div.variables { + let mut c1 = self.clone(); + let mut c2 = div.clone(); + c1.unify_variables(&mut c2); + return self.quot_rem(&c2, abort_on_remainder); + } + if self.nterms() == div.nterms() { if self == div { return (self.one(), self.zero()); } // check if one is a multiple of the other - let (q, r) = self.field.quot_rem(&self.lcoeff(), &div.lcoeff()); + let (q, r) = self.ring.quot_rem(&self.lcoeff(), &div.lcoeff()); if F::is_zero(&r) && self @@ -2512,7 +2567,7 @@ impl MultivariatePolynomial { && self .into_iter() .zip(div) - .all(|(t1, t2)| &self.field.mul(t2.coefficient, &q) == t1.coefficient) + .all(|(t1, t2)| &self.ring.mul(t2.coefficient, &q) == t1.coefficient) { return (self.constant(q), self.zero()); } @@ -2536,7 +2591,7 @@ impl MultivariatePolynomial { } for c in &mut q.coefficients { - let (quot, rem) = q.field.quot_rem(c, dive.coefficient); + let (quot, rem) = q.ring.quot_rem(c, dive.coefficient); *c = quot; if !F::is_zero(&rem) { // TODO: support upgrade to a RationalField @@ -2552,7 +2607,7 @@ impl MultivariatePolynomial { .map(|i| self.degree(i).to_u32() as usize + div.degree(i).to_u32() as usize) .collect(); - if div.field.is_one(&div.lcoeff()) && degree_sum.iter().filter(|x| **x > 0).count() == 1 { + if div.ring.is_one(&div.lcoeff()) && degree_sum.iter().filter(|x| **x > 0).count() == 1 { return self.quot_rem_univariate_monic(div); } @@ -2576,7 +2631,7 @@ impl MultivariatePolynomial { /// Heap division for multivariate polynomials, using a cache so that only unique /// monomial exponents appear in the heap. /// Reference: "Sparse polynomial division using a heap" by Monagan, Pearce (2011) - pub fn heap_division( + fn heap_division( &self, div: &MultivariatePolynomial, abort_on_remainder: bool, @@ -2612,7 +2667,7 @@ impl MultivariatePolynomial { for (s, e) in m.iter_mut().zip(h.peek().unwrap().as_slice()) { *s = *e; } - c = self.field.zero(); + c = self.ring.zero(); } if let Some(monomial) = h.peek() { @@ -2621,7 +2676,7 @@ impl MultivariatePolynomial { let mut qs = cache.remove(&m).unwrap(); for (i, j, next_in_divisor) in qs.drain(..) { - self.field.sub_mul_assign( + self.ring.sub_mul_assign( &mut c, &q.coefficients[i], div.coefficient_back(j), @@ -2715,7 +2770,7 @@ impl MultivariatePolynomial { } if div.last_exponents().iter().zip(&m).all(|(ge, me)| me >= ge) { - let (quot, rem) = self.field.quot_rem(&c, &div.lcoeff()); + let (quot, rem) = self.ring.quot_rem(&c, &div.lcoeff()); if !F::is_zero(&rem) { if abort_on_remainder { r = self.one(); @@ -2884,7 +2939,7 @@ impl MultivariatePolynomial { k += 1; } else { m = *h.peek().unwrap(); - c = self.field.zero(); + c = self.ring.zero(); } if let Some(monomial) = h.peek() { @@ -2894,7 +2949,7 @@ impl MultivariatePolynomial { let mut qs = cache.remove(&m).unwrap(); for (i, j, next_in_divisor) in qs.drain(..) { // TODO: use fraction-free routines - self.field.sub_mul_assign( + self.ring.sub_mul_assign( &mut c, &q.coefficients[i], div.coefficient_back(j), @@ -2971,7 +3026,7 @@ impl MultivariatePolynomial { let q_e = divides(m, pack_div[pack_div.len() - 1], pack_u8); if let Some(q_e) = q_e { - let (quot, rem) = self.field.quot_rem(&c, &div.lcoeff()); + let (quot, rem) = self.ring.quot_rem(&c, &div.lcoeff()); if !F::is_zero(&rem) { if abort_on_remainder { r = self.one(); @@ -3082,6 +3137,13 @@ impl MultivariatePolynomial { /// Compute the p-adic expansion of the polynomial. /// It returns `[a0, a1, a2, ...]` such that `a0 + a1 * p^1 + a2 * p^2 + ... = self`. pub fn p_adic_expansion(&self, p: &Self) -> Vec { + if self.variables != p.variables { + let mut c1 = self.clone(); + let mut c2 = p.clone(); + c1.unify_variables(&mut c2); + return c1.p_adic_expansion(&c2); + } + let mut res = vec![]; let mut r = self.clone(); while !r.is_zero() { @@ -3097,8 +3159,8 @@ impl MultivariatePolynomial { /// Make the polynomial monic, i.e., make the leading coefficient `1` by /// multiplying all monomials with `1/lcoeff`. pub fn make_monic(self) -> Self { - if self.lcoeff() != self.field.one() { - let ci = self.field.inv(&self.lcoeff()); + if self.lcoeff() != self.ring.one() { + let ci = self.ring.inv(&self.lcoeff()); self.mul_coeff(ci) } else { self @@ -3120,10 +3182,7 @@ impl MultivariatePolynomial { exp.copy_from_slice(x.exponents); let pow = exp[var].to_u32() as u64; exp[var] += E::one(); - res.append_monomial( - self.field.div(x.coefficient, &self.field.nth(pow + 1)), - &exp, - ); + res.append_monomial(self.ring.div(x.coefficient, &self.ring.nth(pow + 1)), &exp); } res @@ -3146,12 +3205,12 @@ impl MultivariatePolynomial { if div.nterms() == 1 { // calculate inverse once - let inv = self.field.inv(&div.coefficients[0]); + let inv = self.ring.inv(&div.coefficients[0]); if div.is_constant() { let mut q = self.clone(); for c in &mut q.coefficients { - self.field.mul_assign(c, &inv); + self.ring.mul_assign(c, &inv); } return (q, self.zero()); @@ -3163,7 +3222,7 @@ impl MultivariatePolynomial { for m in self.into_iter() { if m.exponents.iter().zip(dive).all(|(a, b)| a >= b) { - q.coefficients.push(self.field.mul(m.coefficient, &inv)); + q.coefficients.push(self.ring.mul(m.coefficient, &inv)); for (ee, ed) in m.exponents.iter().zip(dive) { q.exponents.push(*ee - *ed); @@ -3177,22 +3236,22 @@ impl MultivariatePolynomial { } // normalize the lcoeff to 1 to prevent a costly inversion - if !self.field.is_one(&div.lcoeff()) { + if !self.ring.is_one(&div.lcoeff()) { let o = div.lcoeff(); - let inv = self.field.inv(&div.lcoeff()); + let inv = self.ring.inv(&div.lcoeff()); for c in &mut div.coefficients { - self.field.mul_assign(c, &inv); + self.ring.mul_assign(c, &inv); } let mut res = self.quot_rem_univariate_monic(div); for c in &mut res.0.coefficients { - self.field.mul_assign(c, &inv); + self.ring.mul_assign(c, &inv); } for c in &mut div.coefficients { - self.field.mul_assign(c, &o); + self.ring.mul_assign(c, &o); } return res; @@ -3228,10 +3287,10 @@ impl MultivariatePolynomial { pub fn eea_univariate(&self, other: &Self) -> (Self, Self, Self) { let mut r0 = self.clone().make_monic(); let mut r1 = other.clone().make_monic(); - let mut s0 = self.constant(self.field.inv(&self.lcoeff())); + let mut s0 = self.constant(self.ring.inv(&self.lcoeff())); let mut s1 = self.zero(); let mut t0 = self.zero(); - let mut t1 = self.constant(self.field.inv(&other.lcoeff())); + let mut t1 = self.constant(self.ring.inv(&other.lcoeff())); while !r1.is_zero() { let (q, r) = r0.quot_rem_univariate(&mut r1); @@ -3239,7 +3298,7 @@ impl MultivariatePolynomial { return (r1, s1, t1); } - let a = self.field.inv(&r.lcoeff()); + let a = self.ring.inv(&r.lcoeff()); (r1, r0) = (r.mul_coeff(a.clone()), r1); (s1, s0) = ((s0 - &q * &s1).mul_coeff(a.clone()), s1); (t1, t0) = ((t0 - q * &t1).mul_coeff(a), t1); @@ -3289,7 +3348,7 @@ impl MultivariatePolynomial { let mut exp = self.last_exponents().to_vec(); exp[var] = E::from_u32(deg_n) + E::from_u32(deg_d) + E::one(); - let mut v0 = self.monomial(self.field.one(), exp); + let mut v0 = self.monomial(self.ring.one(), exp); let mut v1 = self.zero(); let mut w0 = self.clone(); @@ -3315,11 +3374,11 @@ impl MultivariatePolynomial { let y_poly = self.to_univariate_polynomial_list(var); let mut sample_powers = Vec::with_capacity(d + 1); - let mut accum = self.field.one(); + let mut accum = self.ring.one(); - sample_powers.push(self.field.one()); + sample_powers.push(self.ring.one()); for _ in 0..d { - self.field.mul_assign(&mut accum, shift); + self.ring.mul_assign(&mut accum, shift); sample_powers.push(accum.clone()); } let mut v = vec![self.zero(); d + 1]; @@ -3335,8 +3394,8 @@ impl MultivariatePolynomial { } let mut poly = self.zero(); - let mut accum_inv = self.field.one(); - let sample_point_inv = self.field.inv(shift); + let mut accum_inv = self.ring.one(); + let sample_point_inv = self.ring.inv(shift); for (i, mut v) in v.into_iter().enumerate() { v = v.mul_coeff(accum_inv.clone()); @@ -3348,7 +3407,7 @@ impl MultivariatePolynomial { poly.append_monomial(m.coefficient.clone(), m.exponents); } - self.field.mul_assign(&mut accum_inv, &sample_point_inv); + self.ring.mul_assign(&mut accum_inv, &sample_point_inv); } poly @@ -3359,7 +3418,7 @@ impl MultivariatePolynomial, E> { /// Convert the polynomial to a multivariate polynomial that contains the /// variable in the number field. pub fn from_number_field(&self) -> MultivariatePolynomial { - let var = &self.field.poly().get_vars_ref()[0]; + let var = &self.ring.poly().get_vars_ref()[0]; let var_map = if !self.get_vars_ref().contains(var) { let mut v = self.get_vars_ref().to_vec(); @@ -3372,7 +3431,7 @@ 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); + MultivariatePolynomial::new(&self.ring.poly().ring, 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); @@ -3393,7 +3452,7 @@ impl From<&MultivariatePolynomial> MultivariatePolynomial { coefficients: val.coefficients.iter().map(|x| x.into()).collect(), exponents: val.exponents.clone(), - field: Q, + ring: Q, variables: val.variables.clone(), _phantom: PhantomData, } @@ -3448,7 +3507,7 @@ impl<'a, F: Ring, E: Exponent, O: MonomialOrder> IntoIterator #[cfg(test)] mod test { - use crate::{atom::Atom, domains::integer::Z}; + use crate::{atom::Atom, domains::integer::Z, state::State}; #[test] fn mul_packed() { @@ -3518,4 +3577,28 @@ mod test { Atom::parse("1-v8-v8*v9-v7-v6-v5-v4+v3-4*v2+v2*v3^2+v2^2*v3").unwrap() ); } + + #[test] + fn fuse_variables() { + let p1 = Atom::parse("v1+v2") + .unwrap() + .to_polynomial::<_, u8>(&Z, None); + let p2 = Atom::parse("v4").unwrap().to_polynomial::<_, u8>(&Z, None); + + let p3 = Atom::parse("v3") + .unwrap() + .to_polynomial::<_, u8>(&Z, p1.variables.clone().into()); + + let r = p1 * &p2 + p3; + + assert_eq!( + r.get_vars_ref(), + &[ + State::get_symbol("v1").into(), + State::get_symbol("v2").into(), + State::get_symbol("v4").into(), + State::get_symbol("v3").into() + ] + ); + } } diff --git a/src/poly/univariate.rs b/src/poly/univariate.rs index 4b40c0e4..6f8880ed 100644 --- a/src/poly/univariate.rs +++ b/src/poly/univariate.rs @@ -560,7 +560,10 @@ impl Add for UnivariatePolynomial { fn add(mut self, mut other: Self) -> Self::Output { assert_eq!(self.field, other.field); - assert_eq!(self.variable, other.variable); + + if self.variable != other.variable { + panic!("Cannot multiply polynomials with different variables"); + } if self.is_zero() { return other; @@ -623,10 +626,16 @@ impl<'a, 'b, F: Ring> Mul<&'a UnivariatePolynomial> for &'b UnivariatePolynom #[inline] fn mul(self, rhs: &'a UnivariatePolynomial) -> Self::Output { + assert_eq!(self.field, rhs.field); + if self.is_zero() || rhs.is_zero() { return self.zero(); } + if self.variable != rhs.variable { + panic!("Cannot multiply polynomials with different variables"); + } + let n = self.degree(); let m = rhs.degree(); @@ -736,6 +745,10 @@ impl UnivariatePolynomial { return Some(self.clone()); } + if self.variable != div.variable { + panic!("Cannot divide with different variables"); + } + // check if the leading coefficients divide if !F::is_zero(&self.field.rem(&self.lcoeff(), &div.lcoeff())) { return None; @@ -799,6 +812,10 @@ impl UnivariatePolynomial { return (self.clone(), self.clone()); } + if self.variable != div.variable { + panic!("Cannot divide with different variables"); + } + let mut n = self.degree(); let m = div.degree(); @@ -837,6 +854,10 @@ impl UnivariatePolynomial { /// Compute the p-adic expansion of the polynomial. /// It returns `[a0, a1, a2, ...]` such that `a0 + a1 * p^1 + a2 * p^2 + ... = self`. pub fn p_adic_expansion(&self, p: &Self) -> Vec { + if self.variable != p.variable { + panic!("Cannot apply p-adic expansion with different variables"); + } + let mut res = vec![]; let mut r = self.clone(); while !r.is_zero() { @@ -915,6 +936,10 @@ impl UnivariatePolynomial { /// Compute `(g, s, t)` where `self * s + other * t = g` /// by means of the extended Euclidean algorithm. pub fn eea(&self, other: &Self) -> (Self, Self, Self) { + if self.variable != other.variable { + panic!("Cannot apply EEA with different variables"); + } + let mut r0 = self.clone().make_monic(); let mut r1 = other.clone().make_monic(); let mut s0 = self.constant(self.field.inv(&self.lcoeff())); @@ -978,6 +1003,10 @@ impl UnivariatePolynomial { return self.clone(); } + if self.variable != b.variable { + panic!("Cannot compute GCD of polynomials with different variables"); + } + let mut c = self.clone(); let mut d = b.clone(); if self.degree() < b.degree() { @@ -1012,6 +1041,10 @@ impl UnivariatePolynomial { return (self.clone(), self.clone()); } + if self.variable != div.variable { + panic!("Cannot divide polynomials with different variables"); + } + let mut n = self.degree(); let m = div.degree(); @@ -1038,8 +1071,11 @@ impl UnivariatePolynomial { impl UnivariatePolynomial> { /// Convert a univariate polynomial of multivariate polynomials to a multivariate polynomial. pub fn flatten(self) -> MultivariatePolynomial { - let Some(pos) = self - .field + if self.is_zero() { + return self.field.zero(); + } + + let Some(pos) = self.coefficients[0] .variables .iter() .position(|x| x == self.variable.as_ref()) @@ -1047,18 +1083,15 @@ impl UnivariatePolynomial> { panic!("Variable not found in the field"); }; + let n_vars = self.coefficients[0].get_vars().len(); let mut res = MultivariatePolynomial::new( &self.field.ring, self.degree().into(), - self.field.variables.clone(), + self.coefficients[0].get_vars().clone(), ); for (p, mut c) in self.coefficients.into_iter().enumerate() { - for (e, nc) in c - .exponents - .chunks_mut(self.field.variables.len()) - .zip(c.coefficients) - { + for (e, nc) in c.exponents.chunks_mut(n_vars).zip(c.coefficients) { e[pos] = E::from_u32(p as u32); res.append_monomial(nc, e); } diff --git a/src/printer.rs b/src/printer.rs index cc7e7edd..b541e42f 100644 --- a/src/printer.rs +++ b/src/printer.rs @@ -870,12 +870,12 @@ impl<'a, R: Ring, E: Exponent> Display for FactorizedRationalPolynomialPrinter<' fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result { if self.opts.explicit_rational_polynomial { if !R::is_zero(&self.poly.numer_coeff) - && !self.poly.numerator.field.is_one(&self.poly.numer_coeff) + && !self.poly.numerator.ring.is_one(&self.poly.numer_coeff) { f.write_fmt(format_args!( "[{}]*", RingPrinter { - ring: &self.poly.numerator.field, + ring: &self.poly.numerator.ring, element: &self.poly.numer_coeff, opts: self.opts, in_product: false @@ -884,7 +884,7 @@ impl<'a, R: Ring, E: Exponent> Display for FactorizedRationalPolynomialPrinter<' } if self.poly.denominators.is_empty() - && self.poly.numerator.field.is_one(&self.poly.denom_coeff) + && self.poly.numerator.ring.is_one(&self.poly.denom_coeff) { if self.poly.numerator.is_zero() { f.write_char('0')?; @@ -907,11 +907,11 @@ impl<'a, R: Ring, E: Exponent> Display for FactorizedRationalPolynomialPrinter<' }, ))?; - if !self.poly.numerator.field.is_one(&self.poly.denom_coeff) { + if !self.poly.numerator.ring.is_one(&self.poly.denom_coeff) { f.write_fmt(format_args!( ",{},1", RingPrinter { - ring: &self.poly.numerator.field, + ring: &self.poly.numerator.ring, element: &self.poly.denom_coeff, opts: self.opts, in_product: false @@ -942,13 +942,13 @@ impl<'a, R: Ring, E: Exponent> Display for FactorizedRationalPolynomialPrinter<' } if self.poly.denominators.is_empty() - && self.poly.numerator.field.is_one(&self.poly.denom_coeff) + && self.poly.numerator.ring.is_one(&self.poly.denom_coeff) { - if !self.poly.numerator.field.is_one(&self.poly.numer_coeff) { + if !self.poly.numerator.ring.is_one(&self.poly.numer_coeff) { f.write_fmt(format_args!( "{}", RingPrinter { - ring: &self.poly.numerator.field, + ring: &self.poly.numerator.ring, element: &self.poly.numer_coeff, opts: self.opts, in_product: false @@ -956,10 +956,10 @@ impl<'a, R: Ring, E: Exponent> Display for FactorizedRationalPolynomialPrinter<' ))?; } - if (self.poly.numerator.field.is_one(&self.poly.numer_coeff) && !self.add_parentheses) + if (self.poly.numerator.ring.is_one(&self.poly.numer_coeff) && !self.add_parentheses) || self.poly.numerator.nterms() < 2 { - if !self.poly.numerator.field.is_one(&self.poly.numer_coeff) { + if !self.poly.numerator.ring.is_one(&self.poly.numer_coeff) { if self.poly.numerator.is_one() { return Ok(()); } @@ -976,7 +976,7 @@ impl<'a, R: Ring, E: Exponent> Display for FactorizedRationalPolynomialPrinter<' } )) } else { - if !self.poly.numerator.field.is_one(&self.poly.numer_coeff) { + if !self.poly.numerator.ring.is_one(&self.poly.numer_coeff) { if self.poly.numerator.is_one() { return Ok(()); } @@ -995,11 +995,11 @@ impl<'a, R: Ring, E: Exponent> Display for FactorizedRationalPolynomialPrinter<' } } else { if self.opts.latex { - if !self.poly.numerator.field.is_one(&self.poly.numer_coeff) { + if !self.poly.numerator.ring.is_one(&self.poly.numer_coeff) { f.write_fmt(format_args!( "{} ", RingPrinter { - ring: &self.poly.numerator.field, + ring: &self.poly.numerator.ring, element: &self.poly.numer_coeff, opts: self.opts, in_product: false @@ -1016,11 +1016,11 @@ impl<'a, R: Ring, E: Exponent> Display for FactorizedRationalPolynomialPrinter<' }, ))?; - if !self.poly.numerator.field.is_one(&self.poly.denom_coeff) { + if !self.poly.numerator.ring.is_one(&self.poly.denom_coeff) { f.write_fmt(format_args!( "{}", RingPrinter { - ring: &self.poly.numerator.field, + ring: &self.poly.numerator.ring, element: &self.poly.denom_coeff, opts: self.opts, in_product: false @@ -1052,11 +1052,11 @@ impl<'a, R: Ring, E: Exponent> Display for FactorizedRationalPolynomialPrinter<' return f.write_str("}}"); } - if !self.poly.numerator.field.is_one(&self.poly.numer_coeff) { + if !self.poly.numerator.ring.is_one(&self.poly.numer_coeff) { f.write_fmt(format_args!( "{}*", RingPrinter { - ring: &self.poly.numerator.field, + ring: &self.poly.numerator.ring, element: &self.poly.numer_coeff, opts: self.opts, in_product: false @@ -1088,7 +1088,7 @@ impl<'a, R: Ring, E: Exponent> Display for FactorizedRationalPolynomialPrinter<' return f.write_fmt(format_args!( "{}", RingPrinter { - ring: &self.poly.numerator.field, + ring: &self.poly.numerator.ring, element: &self.poly.denom_coeff, opts: self.opts, in_product: true @@ -1096,7 +1096,7 @@ impl<'a, R: Ring, E: Exponent> Display for FactorizedRationalPolynomialPrinter<' )); } - if self.poly.numerator.field.is_one(&self.poly.denom_coeff) + if self.poly.numerator.ring.is_one(&self.poly.denom_coeff) && self.poly.denominators.len() == 1 && self.poly.denominators[0].0.nterms() == 1 && self.poly.denominators[0].1 == 1 @@ -1104,7 +1104,7 @@ impl<'a, R: Ring, E: Exponent> Display for FactorizedRationalPolynomialPrinter<' let (d, _) = &self.poly.denominators[0]; let var_count = d.exponents.iter().filter(|x| !x.is_zero()).count(); - if var_count == 0 || d.field.is_one(&d.coefficients[0]) && var_count == 1 { + if var_count == 0 || d.ring.is_one(&d.coefficients[0]) && var_count == 1 { return f.write_fmt(format_args!( "{}", PolynomialPrinter { @@ -1117,11 +1117,11 @@ impl<'a, R: Ring, E: Exponent> Display for FactorizedRationalPolynomialPrinter<' f.write_char('(')?; // TODO: add special cases for 1 argument - if !self.poly.numerator.field.is_one(&self.poly.denom_coeff) { + if !self.poly.numerator.ring.is_one(&self.poly.denom_coeff) { f.write_fmt(format_args!( "{}", RingPrinter { - ring: &self.poly.numerator.field, + ring: &self.poly.numerator.ring, element: &self.poly.denom_coeff, opts: self.opts, in_product: true @@ -1283,7 +1283,7 @@ impl<'a, R: Ring, E: Exponent> Display for RationalPolynomialPrinter<'a, R, E> { || self .poly .denominator - .field + .ring .is_one(&self.poly.denominator.coefficients[0]) && var_count == 1 { @@ -1348,26 +1348,26 @@ impl<'a, F: Ring + Display, E: Exponent, O: MonomialOrder> Display let mut is_first_term = true; for monomial in self.poly { let mut is_first_factor = true; - if self.poly.field.is_one(monomial.coefficient) { + if self.poly.ring.is_one(monomial.coefficient) { if !is_first_term { write!(f, "+")?; } } else if monomial .coefficient - .eq(&self.poly.field.neg(&self.poly.field.one())) + .eq(&self.poly.ring.neg(&self.poly.ring.one())) { write!(f, "-")?; } else { if is_first_term { self.poly - .field + .ring .fmt_display(monomial.coefficient, &self.opts, true, f)?; } else { write!( f, "{:+}", RingPrinter { - ring: &self.poly.field, + ring: &self.poly.ring, element: monomial.coefficient, opts: self.opts, in_product: true @@ -1408,7 +1408,7 @@ impl<'a, F: Ring + Display, E: Exponent, O: MonomialOrder> Display } if self.opts.print_finite_field { - Display::fmt(&self.poly.field, f)?; + Display::fmt(&self.poly.ring, f)?; } Ok(()) diff --git a/src/solve.rs b/src/solve.rs index 7ba88fd6..705eeeae 100644 --- a/src/solve.rs +++ b/src/solve.rs @@ -66,7 +66,7 @@ impl<'a> AtomView<'a> { } } - let field = RationalPolynomialField::new(Z, rhs[0].numerator.get_vars()); + let field = RationalPolynomialField::new(Z); let nrows = (mat.len() / rhs.len()) as u32; let m = Matrix::from_linear(mat, nrows, rhs.len() as u32, field.clone()).unwrap(); @@ -163,7 +163,7 @@ mod test { }) .collect(); - let field = RationalPolynomialField::new_from_poly(&rhs_rat[0].numerator); + let field = RationalPolynomialField::from_poly(&rhs_rat[0].numerator); let m = Matrix::from_linear( system_rat, system.len() as u32,