Merge branch 'main' into zero-knowledge

src/field/crandall_field.rs

@@ -11,6 +11,90 @@ use std::iter::{Product, Sum};
 /// EPSILON = 9 * 2**28 - 1
 const EPSILON: u64 = 2415919103;
 
+/// A precomputed 8*8 Cauchy matrix, generated with `Field::mds_8`.
+const CAUCHY_MDS_8: [[CrandallField; 8]; 8] = [
+    [
+        CrandallField(16140901062381928449),
+        CrandallField(2635249153041947502),
+        CrandallField(3074457345215605419),
+        CrandallField(11068046442776179508),
+        CrandallField(13835058053470224385),
+        CrandallField(6148914690431210838),
+        CrandallField(9223372035646816257),
+        CrandallField(1),
+    ],
+    [
+        CrandallField(2049638230143736946),
+        CrandallField(16140901062381928449),
+        CrandallField(2635249153041947502),
+        CrandallField(3074457345215605419),
+        CrandallField(11068046442776179508),
+        CrandallField(13835058053470224385),
+        CrandallField(6148914690431210838),
+        CrandallField(9223372035646816257),
+    ],
+    [
+        CrandallField(5534023221388089754),
+        CrandallField(2049638230143736946),
+        CrandallField(16140901062381928449),
+        CrandallField(2635249153041947502),
+        CrandallField(3074457345215605419),
+        CrandallField(11068046442776179508),
+        CrandallField(13835058053470224385),
+        CrandallField(6148914690431210838),
+    ],
+    [
+        CrandallField(16769767337539665921),
+        CrandallField(5534023221388089754),
+        CrandallField(2049638230143736946),
+        CrandallField(16140901062381928449),
+        CrandallField(2635249153041947502),
+        CrandallField(3074457345215605419),
+        CrandallField(11068046442776179508),
+        CrandallField(13835058053470224385),
+    ],
+    [
+        CrandallField(10760600708254618966),
+        CrandallField(16769767337539665921),
+        CrandallField(5534023221388089754),
+        CrandallField(2049638230143736946),
+        CrandallField(16140901062381928449),
+        CrandallField(2635249153041947502),
+        CrandallField(3074457345215605419),
+        CrandallField(11068046442776179508),
+    ],
+    [
+        CrandallField(5675921252705733081),
+        CrandallField(10760600708254618966),
+        CrandallField(16769767337539665921),
+        CrandallField(5534023221388089754),
+        CrandallField(2049638230143736946),
+        CrandallField(16140901062381928449),
+        CrandallField(2635249153041947502),
+        CrandallField(3074457345215605419),
+    ],
+    [
+        CrandallField(1317624576520973751),
+        CrandallField(5675921252705733081),
+        CrandallField(10760600708254618966),
+        CrandallField(16769767337539665921),
+        CrandallField(5534023221388089754),
+        CrandallField(2049638230143736946),
+        CrandallField(16140901062381928449),
+        CrandallField(2635249153041947502),
+    ],
+    [
+        CrandallField(15987178195121148178),
+        CrandallField(1317624576520973751),
+        CrandallField(5675921252705733081),
+        CrandallField(10760600708254618966),
+        CrandallField(16769767337539665921),
+        CrandallField(5534023221388089754),
+        CrandallField(2049638230143736946),
+        CrandallField(16140901062381928449),
+    ],
+];
+
 /// A field designed for use with the Crandall reduction algorithm.
 ///
 /// Its order is
@@ -224,6 +308,17 @@ impl Field for CrandallField {
         let x74 = x73 * x39;
         x74
     }
+
+    fn mds_8(vec: [Self; 8]) -> [Self; 8] {
+        let mut result = [Self::ZERO; 8];
+        for r in 0..8 {
+            for c in 0..8 {
+                let entry = CAUCHY_MDS_8[r][c];
+                result[r] += entry * vec[c];
+            }
+        }
+        result
+    }
 }
 
 impl Neg for CrandallField {
@@ -268,7 +363,7 @@ impl Sub for CrandallField {
     #[inline]
     #[allow(clippy::suspicious_arithmetic_impl)]
     fn sub(self, rhs: Self) -> Self {
-        let (diff, under) = self.0.overflowing_sub(rhs.0);
+        let (diff, under) = self.0.overflowing_sub(rhs.to_canonical_u64());
         Self(diff.overflowing_add((under as u64) * Self::ORDER).0)
     }
 }

src/field/field.rs

@@ -1,3 +1,4 @@
+use std::convert::TryInto;
 use std::fmt::{Debug, Display};
 use std::hash::Hash;
 use std::iter::{Product, Sum};
@@ -232,6 +233,35 @@ pub trait Field:
         }
     }
 
+    /// Apply an MDS matrix to the given vector. Any MDS matrix can be used, as long as the same one
+    /// is used among calls with the same vector length.
+    ///
+    /// Note that the default implementation is quite slow. If speed is important, this should be
+    /// overridden with a field-specific implementation which applies a precomputed MDS matrix.
+    fn mds(vec: Vec<Self>) -> Vec<Self> {
+        // We use a Cauchy matrix with x_r = n + r, y_c = c.
+        let n = vec.len();
+        let mut result = Vec::with_capacity(n);
+        for r in 0..n {
+            let mut sum = Self::ZERO;
+            for c in 0..n {
+                let x = Self::from_canonical_usize(n + r);
+                let y = Self::from_canonical_usize(c);
+                // This is the (r, c) entry of the Cauchy matrix.
+                let entry = (x - y).inverse();
+                sum += entry * vec[c];
+            }
+            result.push(sum);
+        }
+        result
+    }
+
+    /// Like `mds`, but specialized to n=8. For specific fields, this can be overridden with an
+    /// impl which applies a fast, precomputed 8x8 MDS matrix.
+    fn mds_8(vec: [Self; 8]) -> [Self; 8] {
+        Self::mds(vec.to_vec()).try_into().unwrap()
+    }
+
     fn rand_from_rng<R: Rng>(rng: &mut R) -> Self {
         Self::from_canonical_u64(rng.gen_range(0, Self::ORDER))
     }

src/field/field_testing.rs

@@ -305,6 +305,16 @@ macro_rules! test_arithmetic {
                     }
                 }
             }
+
+            #[test]
+            fn subtraction() {
+                type F = $field;
+
+                let (a, b) = (F::from_canonical_u64((F::ORDER + 1) / 2), F::TWO);
+                let x = a * b;
+                assert_eq!(x, F::ONE);
+                assert_eq!(F::ZERO - x, F::NEG_ONE);
+            }
         }
     };
 }