Add methods for MDS matrices

The defaults are quite slow, but we will override them with fast, precomputed, field-specific matrices; see #4.

src/field/field.rs

@@ -1,11 +1,12 @@
+use std::convert::TryInto;
 use std::fmt::{Debug, Display};
 use std::hash::Hash;
 use std::iter::{Product, Sum};
 use std::ops::{Add, AddAssign, Div, DivAssign, Mul, MulAssign, Neg, Sub, SubAssign};
 
 use num::Integer;
-use rand::rngs::OsRng;
 use rand::Rng;
+use rand::rngs::OsRng;
 
 use crate::util::bits_u64;
 
@@ -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))
     }