Allow for customizable zero test

- Do not do statistical zero test during series expansion

src/derivative.rs

@@ -336,9 +336,16 @@ impl<'a> AtomView<'a> {
             depth.clone()
         };
 
+        // do not do an expensive statistical zero check during the series expansion
+        // TODO: do such a check on the result of the series expansion?
+        let field = AtomField {
+            statistical_zero_test: false,
+            ..Default::default()
+        };
+
         loop {
             let info = Series::new(
-                &AtomField::new(),
+                &field,
                 None,
                 Arc::new(Variable::Symbol(x)),
                 expansion_point.to_owned(),
@@ -472,7 +479,7 @@ impl<'a> AtomView<'a> {
                             let mut current_depth = info.relative_order();
                             while base_series.is_zero() && current_depth < 1000.into() {
                                 let info = Series::new(
-                                    &AtomField::new(),
+                                    info.get_field(),
                                     None,
                                     info.get_variable().clone(),
                                     info.get_expansion_point().clone(),
@@ -557,7 +564,7 @@ impl Mul<&Atom> for &Series<AtomField> {
 
         loop {
             let info = Series::new(
-                &AtomField::new(),
+                self.get_field(),
                 None,
                 self.get_variable().clone(),
                 expansion_point.to_owned(),
@@ -619,7 +626,7 @@ impl Add<&Atom> for &Series<AtomField> {
 
         loop {
             let info = Series::new(
-                &AtomField::new(),
+                self.get_field(),
                 None,
                 self.get_variable().clone(),
                 expansion_point.to_owned(),
@@ -681,7 +688,7 @@ impl Div<&Atom> for &Series<AtomField> {
 
         loop {
             let info = Series::new(
-                &AtomField::new(),
+                self.get_field(),
                 None,
                 self.get_variable().clone(),
                 expansion_point.to_owned(),

src/domains.rs

@@ -143,7 +143,7 @@ pub trait Ring: Clone + PartialEq + Eq + Hash + Debug + Display {
     /// Return the nth element by computing `n * 1`.
     fn nth(&self, n: Integer) -> Self::Element;
     fn pow(&self, b: &Self::Element, e: u64) -> Self::Element;
-    fn is_zero(a: &Self::Element) -> bool;
+    fn is_zero(&self, a: &Self::Element) -> bool;
     fn is_one(&self, a: &Self::Element) -> bool;
     /// Should return `true` iff `gcd(1,x)` returns `1` for any `x`.
     fn one_is_gcd_unit() -> bool;
@@ -317,8 +317,8 @@ impl<R: Ring, C: Clone + Borrow<R>> Ring for WrappedRingElement<R, C> {
         }
     }
 
-    fn is_zero(a: &Self::Element) -> bool {
-        R::is_zero(&a.element)
+    fn is_zero(&self, a: &Self::Element) -> bool {
+        self.ring().is_zero(&a.element)
     }
 
     fn is_one(&self, a: &Self::Element) -> bool {
@@ -410,7 +410,7 @@ impl<R: Ring, C: Clone + Borrow<R>> InternalOrdering for WrappedRingElement<R, C
 
 impl<R: Ring, C: Clone + Borrow<R>> SelfRing for WrappedRingElement<R, C> {
     fn is_zero(&self) -> bool {
-        R::is_zero(&self.element)
+        self.ring().is_zero(&self.element)
     }
 
     fn is_one(&self) -> bool {

src/domains/algebraic_number.rs

@@ -511,7 +511,7 @@ impl<R: EuclideanDomain> Ring for AlgebraicExtension<R> {
         result
     }
 
-    fn is_zero(a: &Self::Element) -> bool {
+    fn is_zero(&self, a: &Self::Element) -> bool {
         a.poly.is_zero()
     }
 

src/domains/atom.rs

@@ -40,6 +40,8 @@ impl<T: Clone + Send + Sync + Fn(AtomView<'_>, &mut Atom) -> bool> Map for T {}
 /// ```
 #[derive(Clone)]
 pub struct AtomField {
+    /// Use a statistical test to check for zero, instead of only checking for an exact zero.
+    pub statistical_zero_test: bool,
     /// Perform a cancellation check of numerators and denominators after a division.
     pub cancel_check_on_division: bool,
     /// A custom normalization function applied after every operation.
@@ -67,6 +69,7 @@ impl Default for AtomField {
 impl AtomField {
     pub fn new() -> AtomField {
         AtomField {
+            statistical_zero_test: true,
             custom_normalization: None,
             cancel_check_on_division: false,
         }
@@ -172,8 +175,12 @@ impl Ring for AtomField {
     }
 
     /// Check if the result could be 0 using a statistical method.
-    fn is_zero(a: &Self::Element) -> bool {
-        !a.as_view().zero_test(10, f64::EPSILON).is_false()
+    fn is_zero(&self, a: &Self::Element) -> bool {
+        if self.statistical_zero_test {
+            !a.as_view().zero_test(10, f64::EPSILON).is_false()
+        } else {
+            a.is_zero()
+        }
     }
 
     fn is_one(&self, a: &Self::Element) -> bool {

src/domains/factorized_rational_polynomial.rs

@@ -166,15 +166,16 @@ impl<R: Ring, E: PositiveExponent> SelfRing for FactorizedRationalPolynomial<R,
         mut state: PrintState,
         f: &mut W,
     ) -> Result<bool, Error> {
-        let has_numer_coeff = !self.numerator.ring.is_one(&self.numer_coeff);
-        let has_denom_coeff = !self.numerator.ring.is_one(&self.denom_coeff);
+        let ring = &self.numerator.ring;
+        let has_numer_coeff = !ring.is_one(&self.numer_coeff);
+        let has_denom_coeff = !ring.is_one(&self.denom_coeff);
 
         if opts.explicit_rational_polynomial {
             if state.in_sum {
                 f.write_char('+')?;
             }
 
-            if !R::is_zero(&self.numer_coeff) && has_numer_coeff {
+            if !ring.is_zero(&self.numer_coeff) && has_numer_coeff {
                 f.write_char('[')?;
                 self.numerator
                     .ring
@@ -218,7 +219,7 @@ impl<R: Ring, E: PositiveExponent> SelfRing for FactorizedRationalPolynomial<R,
             return Ok(false);
         }
 
-        if R::is_zero(&self.numer_coeff) {
+        if ring.is_zero(&self.numer_coeff) {
             if state.in_sum {
                 f.write_str("+0")?;
             } else {
@@ -834,7 +835,7 @@ where
         poly
     }
 
-    fn is_zero(a: &Self::Element) -> bool {
+    fn is_zero(&self, a: &Self::Element) -> bool {
         a.numerator.is_zero()
     }
 

src/domains/finite_field.rs

@@ -374,7 +374,7 @@ impl Ring for Zp {
     }
 
     #[inline]
-    fn is_zero(a: &Self::Element) -> bool {
+    fn is_zero(&self, a: &Self::Element) -> bool {
         a.0 == 0
     }
 
@@ -396,7 +396,7 @@ impl Ring for Zp {
     }
 
     fn try_div(&self, a: &Self::Element, b: &Self::Element) -> Option<Self::Element> {
-        if Self::is_zero(b) {
+        if self.is_zero(b) {
             None
         } else {
             Some(self.div(a, b))
@@ -681,7 +681,7 @@ impl Ring for Zp64 {
     }
 
     #[inline]
-    fn is_zero(a: &Self::Element) -> bool {
+    fn is_zero(&self, a: &Self::Element) -> bool {
         a.0 == 0
     }
 
@@ -703,7 +703,7 @@ impl Ring for Zp64 {
     }
 
     fn try_div(&self, a: &Self::Element, b: &Self::Element) -> Option<Self::Element> {
-        if Self::is_zero(b) {
+        if self.is_zero(b) {
             None
         } else {
             Some(self.div(a, b))
@@ -960,7 +960,7 @@ impl Ring for FiniteField<Two> {
     }
 
     #[inline]
-    fn is_zero(a: &Self::Element) -> bool {
+    fn is_zero(&self, a: &Self::Element) -> bool {
         *a == 0
     }
 
@@ -1218,7 +1218,7 @@ impl Ring for FiniteField<Mersenne64> {
     }
 
     #[inline]
-    fn is_zero(a: &Self::Element) -> bool {
+    fn is_zero(&self, a: &Self::Element) -> bool {
         *a == 0
     }
 
@@ -1240,7 +1240,7 @@ impl Ring for FiniteField<Mersenne64> {
     }
 
     fn try_div(&self, a: &Self::Element, b: &Self::Element) -> Option<Self::Element> {
-        if Self::is_zero(b) {
+        if self.is_zero(b) {
             None
         } else {
             Some(self.div(a, b))
@@ -1491,7 +1491,7 @@ impl Ring for FiniteField<Integer> {
         b.pow(e).symmetric_mod(&self.p)
     }
 
-    fn is_zero(a: &Self::Element) -> bool {
+    fn is_zero(&self, a: &Self::Element) -> bool {
         a.is_zero()
     }
 

src/domains/float.rs

@@ -122,7 +122,7 @@ impl<T: NumericalFloatLike + SingleFloat + Hash + Eq + InternalOrdering> Ring fo
     }
 
     #[inline(always)]
-    fn is_zero(a: &Self::Element) -> bool {
+    fn is_zero(&self, a: &Self::Element) -> bool {
         a.is_zero()
     }
 

src/domains/integer.rs

@@ -1139,7 +1139,7 @@ impl Ring for IntegerRing {
     }
 
     #[inline]
-    fn is_zero(a: &Self::Element) -> bool {
+    fn is_zero(&self, a: &Self::Element) -> bool {
         match a {
             Integer::Natural(r) => *r == 0,
             Integer::Double(_) => false,
@@ -2236,7 +2236,7 @@ impl Ring for MultiPrecisionIntegerRing {
     }
 
     #[inline]
-    fn is_zero(a: &Self::Element) -> bool {
+    fn is_zero(&self, a: &Self::Element) -> bool {
         a.is_zero()
     }
 

src/domains/rational.rs

@@ -298,8 +298,8 @@ impl<R: EuclideanDomain + FractionNormalization> Ring for FractionField<R> {
         }
     }
 
-    fn is_zero(a: &Self::Element) -> bool {
-        R::is_zero(&a.numerator)
+    fn is_zero(&self, a: &Self::Element) -> bool {
+        self.ring.is_zero(&a.numerator)
     }
 
     fn is_one(&self, a: &Self::Element) -> bool {
@@ -320,7 +320,7 @@ impl<R: EuclideanDomain + FractionNormalization> Ring for FractionField<R> {
     }
 
     fn try_div(&self, a: &Self::Element, b: &Self::Element) -> Option<Self::Element> {
-        if Self::is_zero(b) {
+        if self.is_zero(b) {
             None
         } else {
             Some(self.div(a, b))
@@ -484,7 +484,7 @@ impl<R: EuclideanDomain + FractionNormalization> Field for FractionField<R> {
     }
 
     fn inv(&self, a: &Self::Element) -> Self::Element {
-        if R::is_zero(&a.numerator) {
+        if self.ring.is_zero(&a.numerator) {
             panic!("Division by 0");
         }
 

src/domains/rational_polynomial.rs

@@ -653,7 +653,7 @@ where
         poly
     }
 
-    fn is_zero(a: &Self::Element) -> bool {
+    fn is_zero(&self, a: &Self::Element) -> bool {
         a.numerator.is_zero()
     }
 
@@ -1181,7 +1181,7 @@ where
             prs.push(b.clone().make_monic());
             while !prs.last().unwrap().is_zero() {
                 let r = prs[prs.len() - 2].rem(&prs[prs.len() - 1]);
-                if RationalPolynomialField::is_zero(&r.lcoeff()) {
+                if r.lcoeff().is_zero() {
                     break;
                 }
                 prs.push(r.make_monic());

src/poly/factor.rs

@@ -1150,7 +1150,7 @@ where
                     let r = self
                         .ring
                         .sample(&mut rng, (0, characteristic.to_i64().unwrap_or(i64::MAX)));
-                    if !F::is_zero(&r) {
+                    if !self.ring.is_zero(&r) {
                         exp[var] = E::from_u32(i as u32);
                         random_poly.append_monomial(r, &exp);
                     }
@@ -1366,7 +1366,7 @@ where
 
         let mut d = self.degree(interpolation_var).to_u32();
 
-        let shifted_poly = if !F::is_zero(&sample_point) {
+        let shifted_poly = if !self.ring.is_zero(&sample_point) {
             self.shift_var_cached(interpolation_var, &sample_point)
         } else {
             self.clone()
@@ -1435,7 +1435,7 @@ where
 
         rec_factors.push(rest);
 
-        if !F::is_zero(&sample_point) {
+        if !self.ring.is_zero(&sample_point) {
             for x in &mut rec_factors {
                 // shift the polynomial to y - sample
                 *x = x.shift_var_cached(interpolation_var, &self.ring.neg(&sample_point));