Fix field division by handling division by zero (#969)

* fix every place where div was used

* change expect() in mersenne div

* use unwrap_unchecked on evaluate_zerofier

* fix clippy in benches

* update comments

* replace unwrap_unchecked() for expect since still is it possible to get a denominator equal to zero

* fix conflicts

* handle division by zero error

* refactor code to remove unsafe and compute inverse once

* remove unwrap from round_2

* revert commit

---------

Co-authored-by: Nicole <nicole@Nicoles-Air.fibertel.com.ar>

crypto/src/hash/monolith/mod.rs

@@ -166,7 +166,11 @@ impl<const WIDTH: usize, const NUM_FULL_ROUNDS: usize> MonolithMersenne31<WIDTH,
 
         for (i, x_i) in x.iter().enumerate() {
             for (j, yj) in y.iter().enumerate() {
-                output[i] = F::add(&output[i], &F::div(&to_multiply[j], &F::add(x_i, yj)));
+                output[i] = F::add(
+                    &output[i],
+                    // We are using that x_i + yj != 0 because they are both much smaller than the modulus.
+                    &F::div(&to_multiply[j], &F::add(x_i, yj)).unwrap(),
+                );
             }
         }
 

math/benches/fields/baby_bear.rs

@@ -156,7 +156,7 @@ pub fn babybear_u32_ops_benchmarks(c: &mut Criterion) {
         group.bench_with_input(format!("Division {:?}", &i.len()), &i, |bench, i| {
             bench.iter(|| {
                 for (x, y) in i {
-                    black_box(black_box(x) / black_box(y));
+                    black_box(black_box(x) / black_box(y)).unwrap();
                 }
             });
         });
@@ -214,7 +214,7 @@ pub fn babybear_u64_ops_benchmarks(c: &mut Criterion) {
         group.bench_with_input(format!("Division {:?}", &i.len()), &i, |bench, i| {
             bench.iter(|| {
                 for (x, y) in i {
-                    black_box(black_box(x) / black_box(y));
+                    black_box(black_box(x) / black_box(y)).unwrap();
                 }
             });
         });
@@ -266,7 +266,7 @@ pub fn babybear_u32_extension_ops_benchmarks(c: &mut Criterion) {
         group.bench_with_input(format!("Inverse of Fp4 {:?}", &i.len()), &i, |bench, i| {
             bench.iter(|| {
                 for (x, y) in i {
-                    black_box(black_box(x) / black_box(y));
+                    black_box(black_box(x) / black_box(y)).unwrap();
                 }
             });
         });
@@ -328,7 +328,7 @@ pub fn babybear_u64_extension_ops_benchmarks(c: &mut Criterion) {
         group.bench_with_input(format!("Inverse of Fp4 {:?}", &i.len()), &i, |bench, i| {
             bench.iter(|| {
                 for (x, y) in i {
-                    black_box(black_box(x) / black_box(y));
+                    black_box(black_box(x) / black_box(y)).unwrap();
                 }
             });
         });

math/benches/fields/mersenne31.rs

@@ -183,7 +183,7 @@ pub fn mersenne31_ops_benchmarks(c: &mut Criterion) {
         group.bench_with_input(format!("div {:?}", &i.len()), &i, |bench, i| {
             bench.iter(|| {
                 for (x, y) in i {
-                    black_box(black_box(x) / black_box(y));
+                    black_box(black_box(x) / black_box(y)).unwrap();
                 }
             });
         });

math/benches/fields/mersenne31_montgomery.rs

@@ -154,7 +154,7 @@ pub fn mersenne31_mont_ops_benchmarks(c: &mut Criterion) {
         group.bench_with_input(format!("div {:?}", &i.len()), &i, |bench, i| {
             bench.iter(|| {
                 for (x, y) in i {
-                    black_box(black_box(x) / black_box(y));
+                    black_box(black_box(x) / black_box(y)).unwrap();
                 }
             });
         });

math/benches/fields/stark252.rs

@@ -153,7 +153,7 @@ pub fn starkfield_ops_benchmarks(c: &mut Criterion) {
         group.bench_with_input(format!("div {:?}", &i.len()), &i, |bench, i| {
             bench.iter(|| {
                 for (x, y) in i {
-                    black_box(black_box(x) / black_box(y));
+                    black_box(black_box(x) / black_box(y)).unwrap();
                 }
             });
         });

math/benches/fields/u64_goldilocks.rs

@@ -123,7 +123,7 @@ pub fn u64_goldilocks_ops_benchmarks(c: &mut Criterion) {
         group.bench_with_input(format!("div {:?}", &i.len()), &i, |bench, i| {
             bench.iter(|| {
                 for (x, y) in i {
-                    black_box(black_box(x) / black_box(y));
+                    black_box(black_box(x) / black_box(y)).unwrap();
                 }
             });
         });

math/benches/fields/u64_goldilocks_montgomery.rs

@@ -157,7 +157,7 @@ pub fn u64_goldilocks_montgomery_ops_benchmarks(c: &mut Criterion) {
         group.bench_with_input(format!("div {:?}", &i.len()), &i, |bench, i| {
             bench.iter(|| {
                 for (x, y) in i {
-                    black_box(black_box(x) / black_box(y));
+                    black_box(black_box(x) / black_box(y)).unwrap();
                 }
             });
         });

math/src/elliptic_curve/edwards/point.rs

@@ -134,9 +134,12 @@ impl<E: IsEdwards> IsGroup for EdwardsProjectivePoint<E> {
 
         let num_s2 = &y1y2 - E::a() * &x1x2;
         let den_s2 = &one - &dx1x2y1y2;
-
         // SAFETY: The creation of the result point is safe because the inputs are always points that belong to the curve.
-        let point = Self::new([&num_s1 / &den_s1, &num_s2 / &den_s2, one]);
+        // We are using that den_s1 and den_s2 aren't zero.
+        // See Theorem 3.3 from https://eprint.iacr.org/2007/286.pdf.
+        let x_coord = (&num_s1 / &den_s1).unwrap();
+        let y_coord = (&num_s2 / &den_s2).unwrap();
+        let point = Self::new([x_coord, y_coord, one]);
         point.unwrap()
     }
 

math/src/elliptic_curve/montgomery/point.rs

@@ -142,7 +142,11 @@ impl<E: IsMontgomery> IsGroup for MontgomeryProjectivePoint<E> {
                 let x1_square = &x1 * &x1;
                 let num = &x1_square + &x1_square + x1_square + &x1a + x1a + &one;
                 let den = (&b + &b) * &y1;
-                let div = num / den;
+
+                // We are using that den != 0 because b and y1 aren't zero.
+                // b != 0 because the cofficient b of a montgomery elliptic curve has to be different from zero.
+                // y1 != 0 because if not, it woould be the case from above: x2 = x1 and y2 + y1 = 0.
+                let div = unsafe { (num / den).unwrap_unchecked() };
 
                 let new_x = &div * &div * &b - (&x1 + x2) - a;
                 let new_y = div * (x1 - &new_x) - y1;
@@ -156,7 +160,8 @@ impl<E: IsMontgomery> IsGroup for MontgomeryProjectivePoint<E> {
             } else {
                 let num = &y2 - &y1;
                 let den = &x2 - &x1;
-                let div = num / den;
+
+                let div = unsafe { (num / den).unwrap_unchecked() };
 
                 let new_x = &div * &div * E::b() - (&x1 + &x2) - E::a();
                 let new_y = div * (x1 - &new_x) - y1;

math/src/elliptic_curve/short_weierstrass/curves/bls12_377/field_extension.rs

@@ -69,8 +69,9 @@ impl IsField for Degree2ExtensionField {
     }
 
     /// Returns the division of `a` and `b`
-    fn div(a: &Self::BaseType, b: &Self::BaseType) -> Self::BaseType {
-        <Self as IsField>::mul(a, &Self::inv(b).unwrap())
+    fn div(a: &Self::BaseType, b: &Self::BaseType) -> Result<Self::BaseType, FieldError> {
+        let b_inv = &Self::inv(b).map_err(|_| FieldError::DivisionByZero)?;
+        Ok(<Self as IsField>::mul(a, b_inv))
     }
 
     /// Returns a boolean indicating whether `a` and `b` are equal component wise.
@@ -124,9 +125,11 @@ impl IsSubFieldOf<Degree2ExtensionField> for BLS12377PrimeField {
     fn div(
         a: &Self::BaseType,
         b: &<Degree2ExtensionField as IsField>::BaseType,
-    ) -> <Degree2ExtensionField as IsField>::BaseType {
-        let b_inv = Degree2ExtensionField::inv(b).unwrap();
-        <Self as IsSubFieldOf<Degree2ExtensionField>>::mul(a, &b_inv)
+    ) -> Result<<Degree2ExtensionField as IsField>::BaseType, FieldError> {
+        let b_inv = Degree2ExtensionField::inv(b)?;
+        Ok(<Self as IsSubFieldOf<Degree2ExtensionField>>::mul(
+            a, &b_inv,
+        ))
     }
 
     fn sub(
@@ -378,7 +381,7 @@ mod tests {
         let a = Fp6E::from(3);
         let a_extension = Fp12E::from(3);
         let b = Fp12E::from(2);
-        assert_eq!(a / &b, a_extension / b);
+        assert_eq!((a / &b).unwrap(), (a_extension / b).unwrap());
     }
 
     #[test]

math/src/elliptic_curve/short_weierstrass/curves/bls12_381/field_extension.rs

@@ -71,8 +71,9 @@ impl IsField for Degree2ExtensionField {
     }
 
     /// Returns the division of `a` and `b`
-    fn div(a: &Self::BaseType, b: &Self::BaseType) -> Self::BaseType {
-        <Self as IsField>::mul(a, &Self::inv(b).unwrap())
+    fn div(a: &Self::BaseType, b: &Self::BaseType) -> Result<Self::BaseType, FieldError> {
+        let b_inv = &Self::inv(b).map_err(|_| FieldError::DivisionByZero)?;
+        Ok(<Self as IsField>::mul(a, b_inv))
     }
 
     /// Returns a boolean indicating whether `a` and `b` are equal component wise.
@@ -126,9 +127,11 @@ impl IsSubFieldOf<Degree2ExtensionField> for BLS12381PrimeField {
     fn div(
         a: &Self::BaseType,
         b: &<Degree2ExtensionField as IsField>::BaseType,
-    ) -> <Degree2ExtensionField as IsField>::BaseType {
-        let b_inv = Degree2ExtensionField::inv(b).unwrap();
-        <Self as IsSubFieldOf<Degree2ExtensionField>>::mul(a, &b_inv)
+    ) -> Result<<Degree2ExtensionField as IsField>::BaseType, FieldError> {
+        let b_inv = Degree2ExtensionField::inv(b)?;
+        Ok(<Self as IsSubFieldOf<Degree2ExtensionField>>::mul(
+            a, &b_inv,
+        ))
     }
 
     fn sub(
@@ -429,7 +432,7 @@ mod tests {
         let a = FieldElement::<BLS12381PrimeField>::from(3);
         let a_extension = FieldElement::<Degree2ExtensionField>::from(3);
         let b = FieldElement::<Degree2ExtensionField>::from(2);
-        assert_eq!(a / &b, a_extension / b);
+        assert_eq!((a / &b).unwrap(), (a_extension / b).unwrap());
     }
 
     #[test]

math/src/elliptic_curve/short_weierstrass/curves/bn_254/field_extension.rs

@@ -73,8 +73,9 @@ impl IsField for Degree2ExtensionField {
     }
 
     /// Returns the division of `a` and `b`
-    fn div(a: &Self::BaseType, b: &Self::BaseType) -> Self::BaseType {
-        <Self as IsField>::mul(a, &Self::inv(b).unwrap())
+    fn div(a: &Self::BaseType, b: &Self::BaseType) -> Result<Self::BaseType, FieldError> {
+        let b_inv = &Self::inv(b).map_err(|_| FieldError::DivisionByZero)?;
+        Ok(<Self as IsField>::mul(a, b_inv))
     }
 
     /// Returns a boolean indicating whether `a` and `b` are equal component wise.
@@ -128,9 +129,11 @@ impl IsSubFieldOf<Degree2ExtensionField> for BN254PrimeField {
     fn div(
         a: &Self::BaseType,
         b: &<Degree2ExtensionField as IsField>::BaseType,
-    ) -> <Degree2ExtensionField as IsField>::BaseType {
-        let b_inv = Degree2ExtensionField::inv(b).unwrap();
-        <Self as IsSubFieldOf<Degree2ExtensionField>>::mul(a, &b_inv)
+    ) -> Result<<Degree2ExtensionField as IsField>::BaseType, FieldError> {
+        let b_inv = Degree2ExtensionField::inv(b)?;
+        Ok(<Self as IsSubFieldOf<Degree2ExtensionField>>::mul(
+            a, &b_inv,
+        ))
     }
 
     fn sub(
@@ -397,7 +400,7 @@ mod tests {
         let a = Fp6E::from(3);
         let a_extension = Fp12E::from(3);
         let b = Fp12E::from(2);
-        assert_eq!(a / &b, a_extension / b);
+        assert_eq!((a / &b).unwrap(), (a_extension / b).unwrap());
     }
 
     #[test]

math/src/elliptic_curve/short_weierstrass/point.rs

@@ -347,9 +347,22 @@ where
                     z = ByteConversion::from_bytes_le(&bytes[len * 2..])?;
                 }
 
-                let point =
-                    Self::new([x, y, z]).map_err(|_| DeserializationError::FieldFromBytesError)?;
-                Ok(point)
+                let Ok(z_inv) = z.inv() else {
+                    let point = Self::new([x, y, z])
+                        .map_err(|_| DeserializationError::FieldFromBytesError)?;
+                    return if point.is_neutral_element() {
+                        Ok(point)
+                    } else {
+                        Err(DeserializationError::FieldFromBytesError)
+                    };
+                };
+                let x_affine = &x * &z_inv;
+                let y_affine = &y * &z_inv;
+                if E::defining_equation(&x_affine, &y_affine) == FieldElement::zero() {
+                    Self::new([x, y, z]).map_err(|_| DeserializationError::FieldFromBytesError)
+                } else {
+                    Err(DeserializationError::FieldFromBytesError)
+                }
             }
             PointFormat::Uncompressed => {
                 if bytes.len() % 2 != 0 {

math/src/field/element.rs

@@ -329,12 +329,11 @@ where
     F: IsSubFieldOf<L>,
     L: IsField,
 {
-    type Output = FieldElement<L>;
+    type Output = Result<FieldElement<L>, FieldError>;
 
     fn div(self, rhs: &FieldElement<L>) -> Self::Output {
-        Self::Output {
-            value: <F as IsSubFieldOf<L>>::div(&self.value, &rhs.value),
-        }
+        let value = <F as IsSubFieldOf<L>>::div(&self.value, &rhs.value)?;
+        Ok(FieldElement::<L> { value })
     }
 }
 
@@ -343,7 +342,7 @@ where
     F: IsSubFieldOf<L>,
     L: IsField,
 {
-    type Output = FieldElement<L>;
+    type Output = Result<FieldElement<L>, FieldError>;
 
     fn div(self, rhs: FieldElement<L>) -> Self::Output {
         &self / &rhs
@@ -355,7 +354,7 @@ where
     F: IsSubFieldOf<L>,
     L: IsField,
 {
-    type Output = FieldElement<L>;
+    type Output = Result<FieldElement<L>, FieldError>;
 
     fn div(self, rhs: &FieldElement<L>) -> Self::Output {
         &self / rhs
@@ -367,7 +366,7 @@ where
     F: IsSubFieldOf<L>,
     L: IsField,
 {
-    type Output = FieldElement<L>;
+    type Output = Result<FieldElement<L>, FieldError>;
 
     fn div(self, rhs: FieldElement<L>) -> Self::Output {
         self / &rhs

math/src/field/extensions/cubic.rs

@@ -109,8 +109,12 @@ where
     }
 
     /// Returns the division of `a` and `b`
-    fn div(a: &[FieldElement<F>; 3], b: &[FieldElement<F>; 3]) -> [FieldElement<F>; 3] {
-        <Self as IsField>::mul(a, &Self::inv(b).unwrap())
+    fn div(
+        a: &[FieldElement<F>; 3],
+        b: &[FieldElement<F>; 3],
+    ) -> Result<[FieldElement<F>; 3], FieldError> {
+        let b_inv = &Self::inv(b).map_err(|_| FieldError::DivisionByZero)?;
+        Ok(<Self as IsField>::mul(a, b_inv))
     }
 
     /// Returns a boolean indicating whether `a` and `b` are equal component wise.
@@ -180,9 +184,11 @@ where
     fn div(
         a: &Self::BaseType,
         b: &<CubicExtensionField<F, Q> as IsField>::BaseType,
-    ) -> <CubicExtensionField<F, Q> as IsField>::BaseType {
-        let b_inv = <CubicExtensionField<F, Q> as IsField>::inv(b).unwrap();
-        <Self as IsSubFieldOf<CubicExtensionField<F, Q>>>::mul(a, &b_inv)
+    ) -> Result<<CubicExtensionField<F, Q> as IsField>::BaseType, FieldError> {
+        let b_inv = <CubicExtensionField<F, Q> as IsField>::inv(b)?;
+        Ok(<Self as IsSubFieldOf<CubicExtensionField<F, Q>>>::mul(
+            a, &b_inv,
+        ))
     }
 
     fn sub(
@@ -285,15 +291,15 @@ mod tests {
         let a = FEE::new([FE::new(0), FE::new(3), FE::new(2)]);
         let b = FEE::new([-FE::new(2), FE::new(8), FE::new(5)]);
         let expected_result = FEE::new([FE::new(12), FE::new(6), FE::new(1)]);
-        assert_eq!(a / b, expected_result);
+        assert_eq!((a / b).unwrap(), expected_result);
     }
 
     #[test]
     fn test_div_2() {
         let a = FEE::new([FE::new(12), FE::new(5), FE::new(4)]);
         let b = FEE::new([-FE::new(4), FE::new(2), FE::new(2)]);
         let expected_result = FEE::new([FE::new(3), FE::new(8), FE::new(11)]);
-        assert_eq!(a / b, expected_result);
+        assert_eq!((a / b).unwrap(), expected_result);
     }
 
     #[test]
@@ -379,14 +385,14 @@ mod tests {
         let a = FE::new(2);
         let b = FEE::new([-FE::new(2), FE::new(8), FE::new(5)]);
         let expected_result = FEE::new([FE::new(8), FE::new(4), FE::new(10)]);
-        assert_eq!(a / b, expected_result);
+        assert_eq!((a / b).unwrap(), expected_result);
     }
 
     #[test]
     fn test_div_as_subfield_2() {
         let a = FE::new(4);
         let b = FEE::new([-FE::new(4), FE::new(2), FE::new(2)]);
         let expected_result = FEE::new([FE::new(3), FE::new(6), FE::new(11)]);
-        assert_eq!(a / b, expected_result);
+        assert_eq!((a / b).unwrap(), expected_result);
     }
 }