Lde/lookup grand product

cpp/src/barretenberg/honk/composer/standard_honk_composer.test.cpp

@@ -316,6 +316,7 @@ TEST(StandardHonkComposer, SumcheckRelationCorrectness)
 {
     // Create a composer and a dummy circuit with a few gates
     StandardHonkComposer composer = StandardHonkComposer();
+    static const size_t program_width = StandardHonkComposer::program_width;
     fr a = fr::one();
     // Using the public variable to check that public_input_delta is computed and added to the relation correctly
     uint32_t a_idx = composer.add_public_variable(a);
@@ -348,8 +349,9 @@ TEST(StandardHonkComposer, SumcheckRelationCorrectness)
     };
 
     constexpr size_t num_polynomials = honk::StandardArithmetization::NUM_POLYNOMIALS;
-    // Compute grand product polynomial (now all the necessary polynomials are inside the proving key)
-    polynomial z_perm_poly = prover.compute_grand_product_polynomial(beta, gamma);
+    // Compute grand product polynomial
+    polynomial z_perm_poly = prover_library::compute_permutation_grand_product<program_width>(
+        prover.key, prover.wire_polynomials, beta, gamma);
 
     // Create an array of spans to the underlying polynomials to more easily
     // get the transposition.

cpp/src/barretenberg/honk/proof_system/prover.cpp

@@ -1,6 +1,7 @@
 #include "prover.hpp"
 #include <algorithm>
 #include <cstddef>
+#include "barretenberg/honk/proof_system/prover_library.hpp"
 #include "barretenberg/honk/sumcheck/sumcheck.hpp"
 #include <array>
 #include "barretenberg/honk/sumcheck/polynomials/univariate.hpp" // will go away
@@ -87,114 +88,6 @@ template <typename settings> void Prover<settings>::compute_wire_commitments()
     }
 }
 
-/**
- * @brief Compute the permutation grand product polynomial Z_perm(X)
- * *
- * @detail (This description assumes program_width 3). Z_perm may be defined in terms of its values
- * on X_i = 0,1,...,n-1 as Z_perm[0] = 1 and for i = 1:n-1
- *
- *                  (w_1(j) + β⋅id_1(j) + γ) ⋅ (w_2(j) + β⋅id_2(j) + γ) ⋅ (w_3(j) + β⋅id_3(j) + γ)
- * Z_perm[i] = ∏ --------------------------------------------------------------------------------
- *                  (w_1(j) + β⋅σ_1(j) + γ) ⋅ (w_2(j) + β⋅σ_2(j) + γ) ⋅ (w_3(j) + β⋅σ_3(j) + γ)
- *
- * where ∏ := ∏_{j=0:i-1} and id_i(X) = id(X) + n*(i-1). These evaluations are constructed over the
- * course of four steps. For expositional simplicity, write Z_perm[i] as
- *
- *                A_1(j) ⋅ A_2(j) ⋅ A_3(j)
- * Z_perm[i] = ∏ --------------------------
- *                B_1(j) ⋅ B_2(j) ⋅ B_3(j)
- *
- * Step 1) Compute the 2*program_width length-n polynomials A_i and B_i
- * Step 2) Compute the 2*program_width length-n polynomials ∏ A_i(j) and ∏ B_i(j)
- * Step 3) Compute the two length-n polynomials defined by
- *          numer[i] = ∏ A_1(j)⋅A_2(j)⋅A_3(j), and denom[i] = ∏ B_1(j)⋅B_2(j)⋅B_3(j)
- * Step 4) Compute Z_perm[i+1] = numer[i]/denom[i] (recall: Z_perm[0] = 1)
- *
- * Note: Step (4) utilizes Montgomery batch inversion to replace n-many inversions with
- * one batch inversion (at the expense of more multiplications)
- */
-// TODO(#222)(luke): Parallelize
-template <typename settings> Polynomial Prover<settings>::compute_grand_product_polynomial(Fr beta, Fr gamma)
-{
-    using barretenberg::polynomial_arithmetic::copy_polynomial;
-    static const size_t program_width = settings::program_width;
-
-    // Allocate scratch space for accumulators
-    std::array<Fr*, program_width> numerator_accumulator;
-    std::array<Fr*, program_width> denominator_accumulator;
-    for (size_t i = 0; i < program_width; ++i) {
-        numerator_accumulator[i] = static_cast<Fr*>(aligned_alloc(64, sizeof(Fr) * key->circuit_size));
-        denominator_accumulator[i] = static_cast<Fr*>(aligned_alloc(64, sizeof(Fr) * key->circuit_size));
-    }
-
-    // Populate wire and permutation polynomials
-    std::array<std::span<const Fr>, program_width> wires;
-    std::array<std::span<const Fr>, program_width> sigmas;
-    for (size_t i = 0; i < program_width; ++i) {
-        std::string sigma_id = "sigma_" + std::to_string(i + 1) + "_lagrange";
-        wires[i] = wire_polynomials[i];
-        sigmas[i] = key->polynomial_store.get(sigma_id);
-    }
-
-    // Step (1)
-    // TODO(#222)(kesha): Change the order to engage automatic prefetching and get rid of redundant computation
-    for (size_t i = 0; i < key->circuit_size; ++i) {
-        for (size_t k = 0; k < program_width; ++k) {
-            // Note(luke): this idx could be replaced by proper ID polys if desired
-            Fr idx = k * key->circuit_size + i;
-            numerator_accumulator[k][i] = wires[k][i] + (idx * beta) + gamma;            // w_k(i) + β.(k*n+i) + γ
-            denominator_accumulator[k][i] = wires[k][i] + (sigmas[k][i] * beta) + gamma; // w_k(i) + β.σ_k(i) + γ
-        }
-    }
-
-    // Step (2)
-    for (size_t k = 0; k < program_width; ++k) {
-        for (size_t i = 0; i < key->circuit_size - 1; ++i) {
-            numerator_accumulator[k][i + 1] *= numerator_accumulator[k][i];
-            denominator_accumulator[k][i + 1] *= denominator_accumulator[k][i];
-        }
-    }
-
-    // Step (3)
-    for (size_t i = 0; i < key->circuit_size; ++i) {
-        for (size_t k = 1; k < program_width; ++k) {
-            numerator_accumulator[0][i] *= numerator_accumulator[k][i];
-            denominator_accumulator[0][i] *= denominator_accumulator[k][i];
-        }
-    }
-
-    // Step (4)
-    // Use Montgomery batch inversion to compute z_perm[i+1] = numerator_accumulator[0][i] /
-    // denominator_accumulator[0][i]. At the end of this computation, the quotient numerator_accumulator[0] /
-    // denominator_accumulator[0] is stored in numerator_accumulator[0].
-    Fr* inversion_coefficients = &denominator_accumulator[1][0]; // arbitrary scratch space
-    Fr inversion_accumulator = Fr::one();
-    for (size_t i = 0; i < key->circuit_size; ++i) {
-        inversion_coefficients[i] = numerator_accumulator[0][i] * inversion_accumulator;
-        inversion_accumulator *= denominator_accumulator[0][i];
-    }
-    inversion_accumulator = inversion_accumulator.invert(); // perform single inversion per thread
-    for (size_t i = key->circuit_size - 1; i != std::numeric_limits<size_t>::max(); --i) {
-        numerator_accumulator[0][i] = inversion_accumulator * inversion_coefficients[i];
-        inversion_accumulator *= denominator_accumulator[0][i];
-    }
-
-    // Construct permutation polynomial 'z_perm' in lagrange form as:
-    // z_perm = [0 numerator_accumulator[0][0] numerator_accumulator[0][1] ... numerator_accumulator[0][n-2] 0]
-    Polynomial z_perm(key->circuit_size);
-    // We'll need to shift this polynomial to the left by dividing it by X in gemini, so the the 0-th coefficient should
-    // stay zero
-    copy_polynomial(numerator_accumulator[0], &z_perm[1], key->circuit_size - 1, key->circuit_size - 1);
-
-    // free memory allocated for scratch space
-    for (size_t k = 0; k < program_width; ++k) {
-        aligned_free(numerator_accumulator[k]);
-        aligned_free(denominator_accumulator[k]);
-    }
-
-    return z_perm;
-}
-
 /**
  * - Add circuit size, public input size, and public inputs to transcript
  *
@@ -251,7 +144,8 @@ template <typename settings> void Prover<settings>::execute_grand_product_comput
         .public_input_delta = public_input_delta,
     };
 
-    z_permutation = compute_grand_product_polynomial(beta, gamma);
+    z_permutation =
+        prover_library::compute_permutation_grand_product<settings::program_width>(key, wire_polynomials, beta, gamma);
 
     auto commitment = commitment_key->commit(z_permutation);
 

cpp/src/barretenberg/honk/proof_system/prover.hpp

@@ -23,6 +23,7 @@
 #include <utility>
 #include <string>
 #include "barretenberg/honk/pcs/claim.hpp"
+#include "barretenberg/honk/proof_system/prover_library.hpp"
 
 namespace honk {
 
@@ -46,8 +47,6 @@ template <typename settings> class Prover {
 
     void compute_wire_commitments();
 
-    barretenberg::polynomial compute_grand_product_polynomial(Fr beta, Fr gamma);
-
     void construct_prover_polynomials();
 
     plonk::proof& export_proof();