MVPoly: method to generate random multilinear polynomial

mvpoly/src/lib.rs

@@ -1,3 +1,5 @@
+use std::collections::HashMap;
+
 use ark_ff::PrimeField;
 use rand::RngCore;
 
@@ -59,4 +61,60 @@ pub trait MVPoly<F: PrimeField, const N: usize, const D: usize>:
     /// This is a dummy implementation. A cache can be used for the monomials to
     /// speed up the computation.
     fn eval(&self, x: &[F; N]) -> F;
+
+    /// Returns true if the polynomial is homogeneous (of degree `D`).
+    /// As a reminder, a polynomial is homogeneous if all its monomials have the
+    /// same degree.
+    fn is_homogeneous(&self) -> bool;
+
+    /// Evaluate the polynomial at the vector point `x` and the extra variable
+    /// `u` using its homogeneous form of degree D.
+    fn homogeneous_eval(&self, x: &[F; N], u: F) -> F;
+
+    /// Add the monomial `coeff * x_1^{e_1} * ... * x_N^{e_N}` to the
+    /// polynomial, where `e_i` are the values given by the array `exponents`.
+    ///
+    /// For instance, to add the monomial `3 * x_1^2 * x_2^3` to the polynomial,
+    /// one would call `add_monomial([2, 3], 3)`.
+    fn add_monomial(&mut self, exponents: [usize; N], coeff: F);
+
+    /// Compute the cross-terms as described in [Behind Nova: cross-terms
+    /// computation for high degree
+    /// gates](https://hackmd.io/@dannywillems/Syo5MBq90)
+    ///
+    /// The polynomial must not necessarily be homogeneous. For this reason, the
+    /// values `u1` and `u2` represents the extra variable that is used to make
+    /// the polynomial homogeneous.
+    ///
+    /// The homogeneous degree is supposed to be the one defined by the type of
+    /// the polynomial, i.e. `D`.
+    ///
+    /// The output is a map of `D - 1` values that represents the cross-terms
+    /// for each power of `r`.
+    fn compute_cross_terms(
+        &self,
+        eval1: &[F; N],
+        eval2: &[F; N],
+        u1: F,
+        u2: F,
+    ) -> HashMap<usize, F>;
+
+
+    /// Return true if the multi-variate polynomial is multilinear, i.e. if each
+    /// variable in each monomial is of maximum degree 1.
+    fn is_multilinear(&self) -> bool;
+
+    /// Generate a random multilinear polynomial
+    ///
+    /// # Safety
+    ///
+    /// Marked as unsafe to warn the user to use it with caution and to not
+    /// necessarily rely on it for security/randomness in cryptographic
+    /// protocols. The user is responsible for providing its own secure
+    /// polynomial random generator, if needed.
+    ///
+    /// For now, the function is only used for testing.
+    unsafe fn random_multilinear<RNG: RngCore>(_rng: &mut RNG) -> Self {
+        unimplemented!("TODO: use add_monomial")
+    }
 }

mvpoly/src/monomials.rs

@@ -1,4 +1,5 @@
 use ark_ff::{One, PrimeField, Zero};
+use num_integer::binomial;
 use rand::RngCore;
 use std::{
     collections::HashMap,
@@ -8,7 +9,7 @@ use std::{
 
 use crate::{
     prime,
-    utils::{naive_prime_factors, PrimeNumberGenerator},
+    utils::{compute_indices_nested_loop, naive_prime_factors, PrimeNumberGenerator},
     MVPoly,
 };
 
@@ -351,13 +352,6 @@ impl<const N: usize, const D: usize, F: PrimeField> Sparse<F, N, D> {
             .and_modify(|c| *c = coeff)
             .or_insert(coeff);
     }
-
-    pub fn add_monomial(&mut self, exponents: [usize; N], coeff: F) {
-        self.monomials
-            .entry(exponents)
-            .and_modify(|c| *c += coeff)
-            .or_insert(coeff);
-    }
 }
 
 impl<const N: usize, const D: usize, F: PrimeField> MVPoly<F, N, D> for Sparse<F, N, D> {
@@ -440,6 +434,135 @@ impl<const N: usize, const D: usize, F: PrimeField> MVPoly<F, N, D> for Sparse<F
         // Feel free to change
         prime::Dense::random(rng, max_degree).into()
     }
+
+    fn is_homogeneous(&self) -> bool {
+        self.monomials
+            .iter()
+            .all(|(exponents, _)| exponents.iter().sum::<usize>() == D)
+    }
+
+    // IMPROVEME: powers can be cached
+    fn homogeneous_eval(&self, x: &[F; N], u: F) -> F {
+        self.monomials
+            .iter()
+            .map(|(exponents, coeff)| {
+                let mut term = F::one();
+                for (exp, point) in exponents.iter().zip(x.iter()) {
+                    term *= point.pow([*exp as u64]);
+                }
+                term *= u.pow([D as u64 - exponents.iter().sum::<usize>() as u64]);
+                term * coeff
+            })
+            .sum()
+    }
+
+    fn add_monomial(&mut self, exponents: [usize; N], coeff: F) {
+        self.monomials
+            .entry(exponents)
+            .and_modify(|c| *c += coeff)
+            .or_insert(coeff);
+    }
+
+    fn compute_cross_terms(
+        &self,
+        eval1: &[F; N],
+        eval2: &[F; N],
+        u1: F,
+        u2: F,
+    ) -> HashMap<usize, F> {
+        assert!(
+            D >= 2,
+            "The degree of the polynomial must be greater than 2"
+        );
+        let mut cross_terms_by_powers_of_r: HashMap<usize, F> = HashMap::new();
+        // We iterate over each monomial with their respective coefficient
+        // i.e. we do have something like coeff * x_1^d_1 * x_2^d_2 * ... * x_N^d_N
+        self.monomials.iter().for_each(|(exponents, coeff)| {
+            // "Exponents" contains all powers, even the ones that are 0. We must
+            // get rid of them and keep the index to fetch the correct
+            // evaluation later
+            let non_zero_exponents_with_index: Vec<(usize, &usize)> = exponents
+                .iter()
+                .enumerate()
+                .filter(|(_, &d)| d != 0)
+                .collect();
+            // coeff = 0 should not happen as we suppose we have a sparse polynomial
+            // Therefore, skipping a check
+            let non_zero_exponents: Vec<usize> = non_zero_exponents_with_index
+                .iter()
+                .map(|(_, d)| *d)
+                .copied()
+                .collect::<Vec<usize>>();
+            let monomial_degree = non_zero_exponents.iter().sum::<usize>();
+            let u_degree: usize = D - monomial_degree;
+            // Will be used to compute the nested sums
+            // It returns all the indices i_1, ..., i_k for the sums:
+            // Σ_{i_1 = 0}^{n_1} Σ_{i_2 = 0}^{n_2} ... Σ_{i_k = 0}^{n_k}
+            let indices =
+                compute_indices_nested_loop(non_zero_exponents.iter().map(|d| *d + 1).collect());
+            for i in 0..=u_degree {
+                // Add the binomial from the homogeneisation
+                // i.e (u_degree choose i)
+                let u_binomial_term = binomial(u_degree, i);
+                // Now, we iterate over all the indices i_1, ..., i_k, i.e. we
+                // do over the whole sum, and we populate the map depending on
+                // the power of r
+                indices.iter().for_each(|indices| {
+                    let sum_indices = indices.iter().sum::<usize>() + i;
+                    // power of r is Σ (n_k - i_k)
+                    let power_r: usize = D - sum_indices;
+
+                    // If the sum of the indices is 0 or D, we skip the
+                    // computation as the contribution would go in the
+                    // evaluation of the polynomial at each evaluation
+                    // vectors eval1 and eval2
+                    if sum_indices == 0 || sum_indices == D {
+                        return;
+                    }
+                    // Compute
+                    // (n_1 choose i_1) * (n_2 choose i_2) * ... * (n_k choose i_k)
+                    let binomial_term = indices
+                        .iter()
+                        .zip(non_zero_exponents.iter())
+                        .fold(u_binomial_term, |acc, (i, &d)| acc * binomial(d, *i));
+                    let binomial_term = F::from(binomial_term as u64);
+                    // Compute the product x_k^i_k
+                    // We ignore the power as it comes into account for the
+                    // right evaluation.
+                    // NB: we could merge both loops, but we keep them separate
+                    // for readability
+                    let eval_left = indices
+                        .iter()
+                        .zip(non_zero_exponents_with_index.iter())
+                        .fold(F::one(), |acc, (i, (idx, _d))| {
+                            acc * eval1[*idx].pow([*i as u64])
+                        });
+                    // Compute the product x'_k^(n_k - i_k)
+                    let eval_right = indices
+                        .iter()
+                        .zip(non_zero_exponents_with_index.iter())
+                        .fold(F::one(), |acc, (i, (idx, d))| {
+                            acc * eval2[*idx].pow([(*d - *i) as u64])
+                        });
+                    // u1^i * u2^(u_degree - i)
+                    let u = u1.pow([i as u64]) * u2.pow([(u_degree - i) as u64]);
+                    let res = binomial_term * eval_left * eval_right * u;
+                    let res = *coeff * res;
+                    cross_terms_by_powers_of_r
+                        .entry(power_r)
+                        .and_modify(|e| *e += res)
+                        .or_insert(res);
+                })
+            }
+        });
+        cross_terms_by_powers_of_r
+    }
+
+    fn is_multilinear(&self) -> bool {
+        self.monomials
+            .iter()
+            .all(|(exponents, _)| exponents.iter().all(|&d| d <= 1))
+    }
 }
 
 impl<const N: usize, const D: usize, F: PrimeField> From<prime::Dense<F, N, D>>

mvpoly/src/prime.rs

@@ -329,6 +329,85 @@ impl<F: PrimeField, const N: usize, const D: usize> MVPoly<F, N, D> for Dense<F,
                 }
             })
     }
+
+    fn is_homogeneous(&self) -> bool {
+        let normalized_indices = self.normalized_indices.clone();
+        let mut prime_gen = PrimeNumberGenerator::new();
+        let is_homogeneous = normalized_indices
+            .iter()
+            .zip(self.coeff.clone())
+            .all(|(idx, c)| {
+                let decomposition_of_i = naive_prime_factors(*idx, &mut prime_gen);
+                let monomial_degree = decomposition_of_i.iter().fold(0, |acc, (_, d)| acc + d);
+                monomial_degree == D || c == F::zero()
+            });
+        is_homogeneous
+    }
+
+    fn homogeneous_eval(&self, x: &[F; N], u: F) -> F {
+        let normalized_indices = self.normalized_indices.clone();
+        let mut prime_gen = PrimeNumberGenerator::new();
+        let primes = prime_gen.get_first_nth_primes(N);
+        normalized_indices
+            .iter()
+            .zip(self.coeff.clone())
+            .fold(F::zero(), |acc, (idx, c)| {
+                let decomposition_of_i = naive_prime_factors(*idx, &mut prime_gen);
+                let monomial_degree = decomposition_of_i.iter().fold(0, |acc, (_, d)| acc + d);
+                let u_power = D - monomial_degree;
+                let monomial = decomposition_of_i.iter().fold(F::one(), |acc, (p, d)| {
+                    let inv_p = primes.iter().position(|&x| x == *p).unwrap();
+                    let x_p = x[inv_p].pow([*d as u64]);
+                    acc * x_p
+                });
+                acc + c * monomial * u.pow([u_power as u64])
+            })
+    }
+
+    fn add_monomial(&mut self, exponents: [usize; N], coeff: F) {
+        let mut prime_gen = PrimeNumberGenerator::new();
+        let primes = prime_gen.get_first_nth_primes(N);
+        let normalized_index = exponents
+            .iter()
+            .zip(primes.iter())
+            .fold(1, |acc, (d, p)| acc * p.pow(*d as u32));
+        let inv_idx = self
+            .normalized_indices
+            .iter()
+            .position(|&x| x == normalized_index)
+            .unwrap();
+        self.coeff[inv_idx] += coeff;
+    }
+
+    fn compute_cross_terms(
+        &self,
+        _eval1: &[F; N],
+        _eval2: &[F; N],
+        _u1: F,
+        _u2: F,
+    ) -> HashMap<usize, F> {
+        unimplemented!()
+    }
+
+    fn is_multilinear(&self) -> bool {
+        if self.is_zero() {
+            return true;
+        }
+        let normalized_indices = self.normalized_indices.clone();
+        let mut prime_gen = PrimeNumberGenerator::new();
+        normalized_indices
+            .iter()
+            .zip(self.coeff.iter())
+            .all(|(idx, c)| {
+                if c.is_zero() {
+                    true
+                } else {
+                    let decomposition_of_i = naive_prime_factors(*idx, &mut prime_gen);
+                    // Each prime number/variable should appear at most once
+                    decomposition_of_i.iter().all(|(_p, d)| *d <= 1)
+                }
+            })
+    }
 }
 
 impl<F: PrimeField, const N: usize, const D: usize> Dense<F, N, D> {
@@ -411,21 +490,6 @@ impl<F: PrimeField, const N: usize, const D: usize> Dense<F, N, D> {
         normalized_indices
     }
 
-    /// Returns `true` if the polynomial is homoegenous of degree `d`
-    pub fn is_homogeneous(&self) -> bool {
-        let normalized_indices = self.normalized_indices.clone();
-        let mut prime_gen = PrimeNumberGenerator::new();
-        let is_homogeneous = normalized_indices
-            .iter()
-            .zip(self.coeff.clone())
-            .all(|(idx, c)| {
-                let decomposition_of_i = naive_prime_factors(*idx, &mut prime_gen);
-                let monomial_degree = decomposition_of_i.iter().fold(0, |acc, (_, d)| acc + d);
-                monomial_degree == D || c == F::zero()
-            });
-        is_homogeneous
-    }
-
     pub fn increase_degree<const D_PRIME: usize>(&self) -> Dense<F, N, D_PRIME> {
         assert!(D <= D_PRIME, "The degree of the target polynomial must be greater or equal to the degree of the source polynomial");
         let mut result: Dense<F, N, D_PRIME> = Dense::zero();