Introduce mul_by_inverse_unchecked, and use it

src/fields/fp/mod.rs

@@ -421,6 +421,14 @@ impl<F: PrimeField> AllocatedFp<F> {
         other: &Self,
         should_enforce: &Boolean<F>,
     ) -> Result<(), SynthesisError> {
+        // The high level logic is as follows:
+        // We want to check that self - other != 0. We do this by checking that
+        // (self - other).inverse() exists. In more detail, we check the following:
+        // If `should_enforce == true`, then we set `multiplier = (self - other).inverse()`,
+        // and check that (self - other) * multiplier == 1. (i.e., that the inverse exists)
+        //
+        // If `should_enforce == false`, then we set `multiplier == 0`, and check that
+        // (self - other) * 0 == 0, which is always satisfied.
         let multiplier = Self::new_witness(self.cs.clone(), || {
             if should_enforce.value()? {
                 (self.value.get()? - other.value.get()?).inverse().get()

src/fields/mod.rs

@@ -1,5 +1,5 @@
 use ark_ff::{prelude::*, BitIteratorBE};
-use ark_relations::r1cs::SynthesisError;
+use ark_relations::r1cs::{ConstraintSystemRef, SynthesisError};
 use core::{
     fmt::Debug,
     ops::{Add, AddAssign, Mul, MulAssign, Sub, SubAssign},
@@ -158,8 +158,33 @@ pub trait FieldVar<F: Field, ConstraintF: Field>:
     /// Returns `(self / d)`.
     /// The constraint system will be unsatisfiable when `d = 0`.
     fn mul_by_inverse(&self, d: &Self) -> Result<Self, SynthesisError> {
-        let d_inv = d.inverse()?;
-        Ok(d_inv * self)
+        // Enforce that `d` is not zero.
+        d.enforce_not_equal(&Self::zero())?;
+        self.mul_by_inverse_unchecked(d)
+    }
+
+    /// Returns `(self / d)`.
+    /// The constraint system will be unsatisfiable when `d = 0`.
+    /// This method *does not* check if `d == 0`.
+    fn mul_by_inverse_unchecked(&self, d: &Self) -> Result<Self, SynthesisError> {
+        let cs = self.cs().or(d.cs());
+        match cs {
+            // If we're in the constant case, we just allocate a new constant having value equalling
+            // `self * d.inverse()`.
+            ConstraintSystemRef::None => Self::new_constant(
+                cs,
+                self.value()? * d.value()?.inverse().expect("division by zero"),
+            ),
+            // If not, we allocate `result` as a new witness having value `self * d.inverse()`,
+            // and check that `result * d = self`.
+            _ => {
+                let result = Self::new_witness(ark_relations::ns!(cs, "self  * d_inv"), || {
+                    Ok(self.value()? * &d.value()?.inverse().unwrap_or(F::zero())?)
+                })?;
+                result.mul_equals(d, self)?;
+                Ok(result)
+            }
+        }
     }
 
     /// Computes the frobenius map over `self`.

src/groups/curves/short_weierstrass/non_zero_affine.rs

@@ -55,7 +55,9 @@ where
             // y3 = lambda * (x1 - x3) - y1
             let numerator = y2 - y1;
             let denominator = x2 - x1;
-            let lambda = numerator.mul_by_inverse(&denominator)?;
+            // It's okay to use `unchecked` here, because the precondition of `add_unchecked` is that
+            // self != ±other, which means that `numerator` and `denominator` are both non-zero.
+            let lambda = numerator.mul_by_inverse_unchecked(&denominator)?;
             let x3 = lambda.square()? - x1 - x2;
             let y3 = lambda * &(x1 - &x3) - y1;
             Ok(Self::new(x3, y3))
@@ -80,7 +82,9 @@ where
             // y3 = lambda * (x1 - x3) - y1
             let numerator = x1_sqr.double()? + &x1_sqr + P::COEFF_A;
             let denominator = y1.double()?;
-            let lambda = numerator.mul_by_inverse(&denominator)?;
+            // It's okay to use `unchecked` here, because the precondition of `double` is that
+            // self != zero.
+            let lambda = numerator.mul_by_inverse_unchecked(&denominator)?;
             let x3 = lambda.square()? - x1.double()?;
             let y3 = lambda * &(x1 - &x3) - y1;
             Ok(Self::new(x3, y3))
@@ -96,6 +100,7 @@ where
         if [self].is_constant() || other.is_constant() {
             self.double()?.add_unchecked(other)
         } else {
+            // It's okay to use `unchecked` the precondition is that `self != ±other` (i.e. same logic as in `add_unchecked`)
             let (x1, y1) = (&self.x, &self.y);
             let (x2, y2) = (&other.x, &other.y);
 
@@ -105,12 +110,13 @@ where
             // y3 = lambda * (x1 - x3) - y1
             let numerator = y2 - y1;
             let denominator = x2 - x1;
-            let lambda_1 = numerator.mul_by_inverse(&denominator)?;
+            let lambda_1 = numerator.mul_by_inverse_unchecked(&denominator)?;
 
             let x3 = lambda_1.square()? - x1 - x2;
 
             // Calculate final addition slope:
-            let lambda_2 = (lambda_1 + y1.double()?.mul_by_inverse(&(&x3 - x1))?).negate()?;
+            let lambda_2 =
+                (lambda_1 + y1.double()?.mul_by_inverse_unchecked(&(&x3 - x1))?).negate()?;
 
             let x4 = lambda_2.square()? - x1 - x3;
             let y4 = lambda_2 * &(x1 - &x4) - y1;

src/poly/evaluations/univariate/mod.rs

@@ -175,10 +175,11 @@ impl<F: PrimeField> EvaluationsVar<F> {
         let lhs_numerator = &alpha_to_s - &coset_offset_to_size;
         let lhs_denominator = &coset_offset_to_size * FpVar::constant(F::from(self.domain.size()));
 
-        let lhs = lhs_numerator.mul_by_inverse(&lhs_denominator)?;
+        let lhs = lhs_numerator.mul_by_inverse_unchecked(&lhs_denominator)?;
 
         // `rhs` for coset element `a` is \frac{1}{alpha * a^{-1} - 1} = \frac{1}{alpha * offset^{-1} * a'^{-1} - 1}
-        let alpha_coset_offset_inv = interpolation_point.mul_by_inverse(&self.domain.offset)?;
+        let alpha_coset_offset_inv =
+            interpolation_point.mul_by_inverse_unchecked(&self.domain.offset)?;
 
         // `res` stores the sum of all lagrange polynomials evaluated at alpha
         let mut res = FpVar::<F>::zero();
@@ -190,7 +191,7 @@ impl<F: PrimeField> EvaluationsVar<F> {
             debug_assert_eq!(subgroup_points[i] * subgroup_point_inv, F::one());
             // alpha * offset^{-1} * a'^{-1} - 1
             let lag_donom = &alpha_coset_offset_inv * subgroup_point_inv - F::one();
-            let lag_coeff = lhs.mul_by_inverse(&lag_donom)?;
+            let lag_coeff = lhs.mul_by_inverse_unchecked(&lag_donom)?;
 
             let lag_interpoland = &self.evals[i] * lag_coeff;
             res += lag_interpoland