Initial barycentric Langrange interpolation

Implements barycentric Lagrange interpolation. Uses algorithm (3.1) from the
paper "Polynomial Interpolation: Langrange vs Newton" by Wilhelm Werner to find
the barycentric weights, and then evaluates at `Gf256::zero()` using the second
or "true" form of the barycentric interpolation formula.

I also earlier implemented a variant of this algorithm, Algorithm 2, from "A new
efficient algorithm for polynomial interpolation," which uses less total
operations than Werner's version, however, because it uses a lot more
multiplications or divisions (depending on how you choose to write it), it runs
slower given the running time of subtraction/ addition (equal) vs
multiplication, and especially division in the Gf256 module.

The new algorithm takes n^2 / 2 divisions and n^2 subtractions to calculate the
barycentric weights, and another n divisions, n multiplications, and 2n
additions to evaluate the polynomial*. The old algorithm runs in n^2 - n
divisions, n^2 multiplications, and n^2 subtractions. Without knowing the exact
running time of each of these operations, we can't say for sure, but I think a
good guess would be the new algorithm trends toward about 1/3 running time as n
-> infinity. It's also easy to see theoretically that for small n the original
lagrange algorithm is faster. This is backed up by benchmarks, which showed for
n >= 5, the new algorithm is faster. We can see that this is more or less what
we should expect given the running times in n of these algorithms.

To ensure we always run the faster algorithm, I've kept both versions and only
use the new one when 5 or more points are given.

Previously the tests in the lagrange module were allowed to pass nodes to the
interpolation algorithms with x = 0. Genuine shares will not be evaluated at x =
0, since then they would just be the secret, so:

1. Now nodes in tests start at x = 1 like `scheme::secret_share` deals them out.
2. I have added assert statements to reinforce this fact and guard against
   division by 0 panics.

This meant getting rid of the `evaluate_at_works` test, but
`interpolate_evaluate_at_0_eq_evaluate_at` provides a similar test.

Further work will include the use of barycentric weights in the `interpolate`
function.

A couple more interesting things to note about barycentric weights:

* Barycentric weights can be partially computed if less than threshold
  shares are present. When additional shares come in, computation can resume
  with no penalty to the total runtime.
* They can be determined totally independently from the y values of our points,
  and the x value we want to evaluate for. We only need to know the x values of
  our interpolation points.

src/lagrange.rs

@@ -1,12 +1,28 @@
 use gf256::Gf256;
 use poly::Poly;
 
+// Minimum number of points where it becomes faster to use barycentric
+// Lagrange interpolation. Determined by benchmarking.
+const BARYCENTRIC_THRESHOLD: u8 = 5;
+
 /// Evaluates an interpolated polynomial at `Gf256::zero()` where
-/// the polynomial is determined using Lagrangian interpolation
-/// based on the given `points` in the G(2^8) Galois field.
-pub(crate) fn interpolate_at(points: &[(u8, u8)]) -> u8 {
+/// the polynomial is determined using either standard or barycentric
+/// Lagrange interpolation (depending on number of points, `k`) based
+/// on the given `points` in the G(2^8) Galois field.
+pub(crate) fn interpolate_at(k: u8, points: &[(u8, u8)]) -> u8 {
+    if k < BARYCENTRIC_THRESHOLD {
+        lagrange_interpolate_at(points)
+    } else {
+        barycentric_interpolate_at(k as usize, points)
+    }
+}
+
+/// Evaluates the polynomial at `Gf256::zero()` using standard Langrange
+/// interpolation.
+fn lagrange_interpolate_at(points: &[(u8, u8)]) -> u8 {
     let mut sum = Gf256::zero();
     for (i, &(raw_xi, raw_yi)) in points.iter().enumerate() {
+        assert_ne!(raw_xi, 0, "Invalid share x = 0");
         let xi = Gf256::from_byte(raw_xi);
         let yi = Gf256::from_byte(raw_yi);
         let mut prod = Gf256::one();
@@ -23,6 +39,36 @@ pub(crate) fn interpolate_at(points: &[(u8, u8)]) -> u8 {
     sum.to_byte()
 }
 
+/// Barycentric Lagrange interpolation algorithm from "Polynomial
+/// Interpolation: Langrange vs Newton" by Wilhelm Werner. Evaluates
+/// the polynomial at `Gf256::zero()`.
+fn barycentric_interpolate_at(k: usize, points: &[(u8, u8)]) -> u8 {
+    // Compute the barycentric weights `w`.
+    let mut w = vec![Gf256::zero(); k];
+    w[0] = Gf256::one();
+    let mut x = Vec::with_capacity(k);
+    x.push(Gf256::from_byte(points[0].0));
+    for i in 1..k {
+        x.push(Gf256::from_byte(points[i].0));
+        for j in 0..i {
+            let delta = x[j] - x[i];
+            assert_ne!(delta.poly, 0, "Duplicate shares");
+            w[j] /= delta;
+            w[i] -= w[j];
+        }
+    }
+    // Evaluate the second or "true" form of the barycentric
+    // interpolation formula at `Gf256::zero()`.
+    let (mut num, mut denom) = (Gf256::zero(), Gf256::zero());
+    for i in 0..k {
+        assert_ne!(x[i].poly, 0, "Invalid share x = 0");
+        let diff = w[i] / x[i];
+        num += diff * Gf256::from_byte(points[i].1);
+        denom += diff;
+    }
+    (num / denom).to_byte()
+}
+
 /// Computeds the coefficient of the Lagrange polynomial interpolated
 /// from the given `points`, in the G(2^8) Galois field.
 pub(crate) fn interpolate(points: &[(Gf256, Gf256)]) -> Poly {
@@ -31,6 +77,7 @@ pub(crate) fn interpolate(points: &[(Gf256, Gf256)]) -> Poly {
     let mut poly = vec![Gf256::zero(); len];
 
     for &(x, y) in points {
+        assert_ne!(x.poly, 0, "Invalid share x = 0");
         let mut coeffs = vec![Gf256::zero(); len];
         coeffs[0] = y;
 
@@ -71,24 +118,15 @@ mod tests {
 
     quickcheck! {
 
-        fn evaluate_at_works(ys: Vec<u8>) -> TestResult {
-            if ys.is_empty() || ys.len() > std::u8::MAX as usize {
-                return TestResult::discard();
-            }
-
-            let points = ys.iter().enumerate().map(|(x, y)| (x as u8, *y)).collect::<Vec<_>>();
-            let equals = interpolate_at(points.as_slice()) == ys[0];
-
-            TestResult::from_bool(equals)
-        }
-
-
         fn interpolate_evaluate_at_works(ys: Vec<Gf256>) -> TestResult {
             if ys.is_empty() || ys.len() > std::u8::MAX as usize {
                 return TestResult::discard();
             }
 
-            let points = ys.into_iter().enumerate().map(|(x, y)| (gf256!(x as u8), y)).collect::<Vec<_>>();
+            let points = ys.into_iter()
+                           .zip(1..std::u8::MAX)
+                           .map(|(y, x)| (gf256!(x), y))
+                           .collect::<Vec<_>>();
             let poly = interpolate(&points);
 
             for (x, y) in points {
@@ -105,7 +143,10 @@ mod tests {
                 return TestResult::discard();
             }
 
-            let points = ys.into_iter().enumerate().map(|(x, y)| (x as u8, y)).collect::<Vec<_>>();
+            let points = ys.into_iter()
+                           .zip(1..std::u8::MAX)
+                           .map(|(y, x)| (x, y))
+                           .collect::<Vec<_>>();
 
             let elems = points
                 .iter()
@@ -114,7 +155,8 @@ mod tests {
 
             let poly = interpolate(&elems);
 
-            let equals = poly.evaluate_at(Gf256::zero()).to_byte() == interpolate_at(points.as_slice());
+            let equals = poly.evaluate_at(Gf256::zero()).to_byte()
+                == interpolate_at(points.len() as u8, points.as_slice());
 
             TestResult::from_bool(equals)
         }

src/sss/scheme.rs

@@ -101,7 +101,7 @@ impl SSS {
             for s in shares.iter().take(threshold as usize) {
                 col_in.push((s.id, s.data[byteindex]));
             }
-            secret.push(interpolate_at(&*col_in));
+            secret.push(interpolate_at(threshold, &*col_in));
         }
 
         Ok(secret)