opt: change polynomials to the point form to eliminate elliptic curve operations

Author: just-erray

src/cac/vsss.rs

@@ -57,30 +57,125 @@ impl Secp256k1 {
     }
 }
 
-// we use this for both polynomials over scalars and over projective points
+fn precalculated_factorials_and_inverses(n: usize) -> (Vec<Fr>, Vec<Fr>, Vec<Fr>) {
+    let factorial: Vec<Fr> = std::iter::once(Fr::one())
+        .chain((1..n).scan(Fr::one(), |state, i| {
+            *state *= Fr::from(i as u64);
+            Some(*state)
+        }))
+        .collect();
+
+    // inv_fact[i] = 1 / factorial[i]
+    let inv_factorial: Vec<Fr> = (0..n)
+        .rev()
+        .scan(
+            factorial[n - 1]
+                .inverse()
+                .expect("This is guaranteed to be non-zero"),
+            |cur_state, i| {
+                let ith_value = *cur_state;
+                *cur_state *= Fr::from(i as u64);
+                Some(ith_value)
+            },
+        )
+        .collect::<Vec<_>>()
+        .into_iter()
+        .rev()
+        .collect();
+
+    let inv: Vec<Fr> = (0..n)
+        .map(|i| {
+            if i == 0 {
+                Fr::zero() //This should never be used
+            } else {
+                inv_factorial[i] * factorial[i - 1]
+            }
+        })
+        .collect();
+
+    (factorial, inv_factorial, inv)
+}
+
+// This polynomial is in the (point, value) form instead of coefficient form (points are integers in range [0, degree] converted to Fr's)
+// It's used both for polynomials with scalar and projective point domains
 #[derive(Clone, Debug, Serialize, Deserialize)]
 pub struct Polynomial<T>(Vec<T>);
 
 impl<T> Polynomial<T>
 where
-    for<'a> T: Add<T, Output = T> + Mul<&'a Fr, Output = T> + Clone,
+    for<'a> T: Add<T, Output = T> + Mul<&'a Fr, Output = T> + Clone + Zero,
 {
-    // todo: max x an int
-    fn eval_at(&self, x: Fr) -> T {
-        // Horner's method
-        let mut iter = self.0.iter().rev();
-        let mut acc = iter
-            .next()
-            .expect("polynomial must have at least one coefficient")
-            .clone();
-        for coeff in iter {
-            acc = coeff.clone() + acc * &x;
+    #[allow(dead_code)]
+    // naive lagrange interpolation, used for testing
+    fn eval_at(&self, x: usize) -> T {
+        if x < self.0.len() {
+            return self.0[x].clone();
         }
-        acc
+
+        let x_fr = Fr::from(x as u32);
+        self.0
+            .iter()
+            .enumerate()
+            .fold(T::zero(), |result, (i, y_i)| {
+                let x_i = Fr::from(i as u64);
+                // Compute L_i(x)
+                let (num, denum) = self.0.iter().enumerate().filter(|(j, _)| *j != i).fold(
+                    (Fr::one(), Fr::one()),
+                    |(num, denum), (j, _)| {
+                        let x_j = Fr::from(j as u64);
+                        (num * (x_fr - x_j), denum * (x_i - x_j))
+                    },
+                );
+
+                // calculate li = num / denum = num * denum^{-1}
+                let denum_inv = denum.inverse().expect("x_i - x_j must be nonzero");
+                let li = num * denum_inv;
+
+                result + y_i.clone() * &li
+            })
+    }
+
+    /// evaluates the function at smallest consecutive integer points bigger than the degree
+    /// functions similar to [`lagrange_interpolate_whole_polynomial`],
+    fn eval_at_suffix_points(&self, n_points: usize) -> Vec<T> {
+        let n_known = self.0.len();
+        let n = n_known + n_points;
+        let (factorial, inv_factorial, inv) = precalculated_factorials_and_inverses(n);
+
+        // For x, calculates the multiplication of (x - i) for all i in known_points (known_points = 0..=degree)
+        // returns the inverse of the multiplication result, if its one of the known points
+        let get_coeff = |x: usize| {
+            if x < n_known {
+                //inverse
+                let mut result = inv_factorial[x] * inv_factorial[n_known - 1 - x];
+                if (n_known - x).is_multiple_of(2) {
+                    result *= -Fr::one();
+                }
+                result
+            } else {
+                factorial[x] * inv_factorial[x - n_known]
+            }
+        };
+
+        let lagrange_basis_polynomial_coeffs: Vec<Fr> = (0..n_known).map(&get_coeff).collect();
+
+        (n_known..n_known + n_points)
+            .map(|x| {
+                let mut result = T::zero();
+                let nom_coeff = get_coeff(x);
+                for i in 0..n_known {
+                    let whole_coeff: Fr =
+                        lagrange_basis_polynomial_coeffs[i] * inv[x - i] * nom_coeff;
+                    result = result + self.0[i].clone() * &whole_coeff;
+                }
+                result
+            })
+            .collect()
     }
 }
 
 impl Polynomial<Fr> {
+    // Generate with points with Fr coordinates in the range [0, degree] for efficient commitment
     pub fn rand(mut rand: impl Rng, degree: usize) -> Self {
         Self((0..degree + 1).map(|_| Fr::rand(&mut rand)).collect())
     }
@@ -89,11 +184,17 @@ impl Polynomial<Fr> {
         PolynomialCommits(Polynomial(secp.generator_batch_mul(&self.0)))
     }
 
-    // shares are return with 0-based index. However, we evaluate share i at
-    // x = i+1, since the value at x=0 represents the secret
     pub fn shares(&self, num_shares: usize) -> Vec<(usize, Fr)> {
+        let n_known = self.0.len();
+        let suffix_points = self.eval_at_suffix_points(num_shares - n_known);
         (0..num_shares)
-            .map(|i| (i, self.eval_at(Fr::from((i + 1) as u64))))
+            .map(|i| {
+                if i < n_known {
+                    (i, self.0[i])
+                } else {
+                    (i, suffix_points[i - n_known])
+                }
+            })
             .collect()
     }
 
@@ -116,9 +217,15 @@ pub struct ShareCommits(pub Vec<Projective>);
 
 impl ShareCommits {
     pub fn verify(&self, polynomial_commits: &PolynomialCommits) -> Result<(), String> {
+        let n_known = polynomial_commits.0.0.len();
+        let n_unknown = self.0.len() - n_known;
+        let unknown_points = polynomial_commits.0.eval_at_suffix_points(n_unknown);
         for (i, share_commit) in self.0.iter().enumerate() {
-            let recomputed_share_commit = polynomial_commits.0.eval_at(Fr::from((i + 1) as u64));
-
+            let recomputed_share_commit = if i < n_known {
+                polynomial_commits.0.0[i]
+            } else {
+                unknown_points[i - n_known]
+            };
             if share_commit != &recomputed_share_commit {
                 return Err("Share commit verification failed".to_owned());
             }
@@ -161,40 +268,7 @@ pub fn lagrange_interpolate_whole_polynomial(
     assert!(!known_points.is_empty() || !missing_points.is_empty());
 
     let n = known_points.len() + missing_points.len();
-    let factorial: Vec<Fr> = std::iter::once(Fr::one())
-        .chain((1..n).scan(Fr::one(), |state, i| {
-            *state *= Fr::from(i as u64);
-            Some(*state)
-        }))
-        .collect();
-
-    // inv_fact[i] = 1 / factorial[i]
-    let inv_factorial: Vec<Fr> = (0..n)
-        .rev()
-        .scan(
-            factorial[n - 1]
-                .inverse()
-                .expect("This is guaranteed to be non-zero"),
-            |cur_state, i| {
-                let ith_value = *cur_state;
-                *cur_state *= Fr::from(i as u64);
-                Some(ith_value)
-            },
-        )
-        .collect::<Vec<_>>()
-        .into_iter()
-        .rev()
-        .collect();
-
-    let inv: Vec<Fr> = (0..n)
-        .map(|i| {
-            if i == 0 {
-                Fr::zero() //This should never be used
-            } else {
-                inv_factorial[i] * factorial[i - 1]
-            }
-        })
-        .collect();
+    let (factorial, inv_factorial, inv) = precalculated_factorials_and_inverses(n);
 
     // For x, calculates the multiplication of (x - i) for all i in known_points (known_points = 0..n \ missing_points)
     // returns the inverse of the multiplication result, based on the parameter
@@ -252,18 +326,6 @@ mod tests {
     use rand::{SeedableRng, seq::index::sample};
     use rand_chacha::ChaCha20Rng;
     use std::collections::HashSet;
-    #[test]
-    fn test_polynomial_eval() {
-        let polynomial = Polynomial::<Fr>::rand(rand::thread_rng(), 2);
-
-        match polynomial.0.as_slice() {
-            &[a, b, c] => {
-                let x = Fr::rand(&mut rand::thread_rng());
-                assert_eq!(polynomial.eval_at(x), a + b * x + c * x * x);
-            }
-            _ => unreachable!(),
-        }
-    }
 
     #[test]
     fn test_commit_verification() {
@@ -324,7 +386,7 @@ mod tests {
                 .map(|x| x + 1)
                 .collect::<Vec<_>>();
             let polynomial = Polynomial::rand(seed_rng, n_revealed - 1);
-            let points = polynomial.shares(n_total); //points[i].0 = i
+            let points = polynomial.shares(n_total);
 
             let aux_set: HashSet<_> = hidden_points.iter().copied().collect();
             let known_points: Vec<(usize, Fr)> = points
@@ -336,6 +398,23 @@ mod tests {
 
             for (x, y) in hidden_points.into_iter().zip(answer.into_iter()) {
                 assert_eq!(points[x].1, y);
+                assert_eq!(polynomial.eval_at(x), y);
+            }
+        }
+    }
+
+    #[test]
+    fn test_suffix_interpolation() {
+        for (n_revealed, n_hidden) in vec![(5usize, 2usize), (100, 10), (175, 7)] {
+            // Assumes one of the revealed ones is 0, as it will be in application, includes it in the n_revealed ones
+            let n_total = n_revealed + n_hidden;
+            let seed_rng = ChaCha20Rng::seed_from_u64(42);
+            let polynomial = Polynomial::rand(seed_rng, n_revealed - 1);
+            let points = polynomial.shares(n_total);
+            let answer = polynomial.eval_at_suffix_points(n_hidden);
+            for (x, y) in (n_revealed..n_total).into_iter().zip(answer.into_iter()) {
+                assert_eq!(points[x].1, y);
+                assert_eq!(polynomial.eval_at(x), y);
             }
         }
     }