WIP

mvpoly/src/prime.rs

@@ -460,33 +460,68 @@ impl<F: PrimeField, const N: usize, const D: usize> Dense<F, N, D> {
     // IMPROVEME: Dummy implementation, a cache can be used to save the
     // previously computed multiplications and powers.
     // Maybe using a symbolic form could speed up the computation.
-    pub fn compute_cross_terms(&self, eval1: [F; N], eval2: [F; N], r: F, u1: F, u2: F) -> Vec<F> {
+    pub fn compute_cross_terms(
+        &self,
+        eval1: &[F; N],
+        eval2: &[F; N],
+        r: F,
+        u1: F,
+        u2: F,
+    ) -> Vec<F> {
         assert!(
             D >= 2,
             "The degree of the polynomial must be greater than 2"
         );
         let size = D - 1;
         let mut res: Vec<F> = vec![];
         let mut prime_gen = PrimeNumberGenerator::new();
-        let primes = prime_gen.get_first_nth_primes(N);
         // Initializing to zero
         (0..size).for_each(|_| {
             res.push(F::zero());
         });
         // Computing invididual contribution of each monomial
+        // FIXME: handle constant, prime decompo returns empty list of 1.
         self.coeff.iter().enumerate().for_each(|(i, c)| {
             // If the coefficient is zero, we skip the computation as the contribution
             // is null.
-            if *c == F::zero() {
+            if *c == F::zero() || i == 0 {
+                // FIXME: handle constant, prime decompo returns empty list of 1.
                 return;
             } else {
                 let idx = self.normalized_indices[i];
                 let prime_decomposition = naive_prime_factors(idx, &mut prime_gen);
-                let monomial_degree = prime_decomposition.iter().fold(0, |acc, (_, d)| acc + d);
+                // Fetch individual degree
+                let mut degrees: Vec<usize> = prime_decomposition.iter().map(|(_, d)| *d).collect();
+                // No cross-terms
+                if degrees.len() == 1 && (degrees[0] == 0 || degrees[0] == D) {
+                    return;
+                }
+                let monomial_degree: usize = degrees.iter().sum();
                 // We compute the missing degree to homogeneize the polynomial.
                 // The variable u given as a parameter will be used to
                 // homogeneize the polynomial.
-                let homogeneisation_degree = D - monomial_degree;
+                let nb_term = degrees.iter().map(|d| d + 1).product::<usize>();
+                assert!(
+                    nb_term <= 2_usize.pow(D as u32),
+                    "The number of terms must be less than 2^D"
+                );
+                // Simulate loops through all powers
+                (0..nb_term).for_each(|i| {
+                    let mut div = 1;
+                    // Compute indices for the loop, step i
+                    let indices: Vec<usize> = degrees
+                        .iter()
+                        .map(|n_i| {
+                            let k = (i / div) % (n_i + 1);
+                            div = div * (n_i + 1);
+                            k
+                        })
+                        .collect();
+                    assert!(
+                        div == nb_term,
+                        "The division must be equal to the number of terms at the end"
+                    );
+                })
             }
         });
         res

mvpoly/tests/prime.rs

@@ -802,3 +802,15 @@ fn test_mvpoly_mul_pbt() {
     let p2 = unsafe { Dense::<Fp, 4, 6>::random(&mut rng, Some(max_degree)) };
     assert_eq!(p1.clone() * p2.clone(), p2.clone() * p1.clone());
 }
+
+#[test]
+fn test_mvpoly_compute_cross_terms() {
+    let mut rng = o1_utils::tests::make_test_rng(None);
+    let p1 = unsafe { Dense::<Fp, 4, 2>::random(&mut rng, None) };
+    let random_eval1: [Fp; 4] = std::array::from_fn(|_| Fp::rand(&mut rng));
+    let random_eval2: [Fp; 4] = std::array::from_fn(|_| Fp::rand(&mut rng));
+    let u1 = Fp::rand(&mut rng);
+    let u2 = Fp::rand(&mut rng);
+    let r = Fp::rand(&mut rng);
+    p1.compute_cross_terms(&random_eval1, &random_eval2, u1, u2, r);
+}