diff --git a/cpp/src/barretenberg/honk/composer/standard_honk_composer.test.cpp b/cpp/src/barretenberg/honk/composer/standard_honk_composer.test.cpp
index c8e46024be..f45f5279bf 100644
--- a/cpp/src/barretenberg/honk/composer/standard_honk_composer.test.cpp
+++ b/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.
diff --git a/cpp/src/barretenberg/honk/proof_system/prover.cpp b/cpp/src/barretenberg/honk/proof_system/prover.cpp
index 0673260973..9194f054e3 100644
--- a/cpp/src/barretenberg/honk/proof_system/prover.cpp
+++ b/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);
 
diff --git a/cpp/src/barretenberg/honk/proof_system/prover.hpp b/cpp/src/barretenberg/honk/proof_system/prover.hpp
index 20a61cc302..5c9aa39ad8 100644
--- a/cpp/src/barretenberg/honk/proof_system/prover.hpp
+++ b/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();
diff --git a/cpp/src/barretenberg/honk/proof_system/prover.test.cpp b/cpp/src/barretenberg/honk/proof_system/prover.test.cpp
deleted file mode 100644
index 7d615dc26d..0000000000
--- a/cpp/src/barretenberg/honk/proof_system/prover.test.cpp
+++ /dev/null
@@ -1,160 +0,0 @@
-#include "prover.hpp"
-#include "barretenberg/honk/composer/standard_honk_composer.hpp"
-#include "barretenberg/polynomials/polynomial.hpp"
-
-#include "barretenberg/srs/reference_string/file_reference_string.hpp"
-#include <array>
-#include <vector>
-#include <cstddef>
-#include <gtest/gtest.h>
-
-using namespace honk;
-namespace honk_prover_tests {
-
-// field is named Fscalar here because of clash with the Fr
-template <class Fscalar> class ProverTests : public testing::Test {
-
-  public:
-    /**
-     * @brief Test the correctness of the computation of the permutation grand product polynomial z_permutation
-     * @details This test compares a simple, unoptimized, easily readable calculation of the grand product z_permutation
-     * to the optimized implementation used by the prover. It's purpose is to provide confidence that some optimization
-     * introduced into the calculation has not changed the result.
-     */
-    static void test_grand_product_construction()
-    {
-        using barretenberg::polynomial;
-
-        // Define some mock inputs for proving key constructor
-        static const size_t num_gates = 8;
-        static const size_t num_public_inputs = 0;
-        auto reference_string = std::make_shared<bonk::FileReferenceString>(num_gates + 1, "../srs_db/ignition");
-
-        // Instatiate a proving_key and make a pointer to it. This will be used to instantiate a Prover.
-        auto proving_key = std::make_shared<bonk::proving_key>(
-            num_gates, num_public_inputs, reference_string, plonk::ComposerType::STANDARD_HONK);
-
-        static const size_t program_width = StandardProver::settings_::program_width;
-
-        // Construct mock wire and permutation polynomials and add them to the proving_key.
-        // Note: for the purpose of checking the consistency between two methods of computing z_perm, these polynomials
-        // can simply be random. We're not interested in the particular properties of the result.
-        std::vector<polynomial> wires;
-        std::vector<polynomial> sigmas;
-        for (size_t i = 0; i < program_width; ++i) {
-            polynomial wire_poly(proving_key->circuit_size);
-            polynomial sigma_poly(proving_key->circuit_size);
-            for (size_t j = 0; j < proving_key->circuit_size; ++j) {
-                wire_poly[j] = Fscalar::random_element();
-                sigma_poly[j] = Fscalar::random_element();
-            }
-            // Copy the wires and sigmas for use in constructing the test-owned z_perm
-            wires.emplace_back(wire_poly);
-            sigmas.emplace_back(sigma_poly);
-
-            // Add polys to proving_key; to be used by the prover in constructing it's own z_perm
-            std::string wire_id = "w_" + std::to_string(i + 1) + "_lagrange";
-            std::string sigma_id = "sigma_" + std::to_string(i + 1) + "_lagrange";
-            proving_key->polynomial_store.put(wire_id, std::move(wire_poly));
-            proving_key->polynomial_store.put(sigma_id, std::move(sigma_poly));
-        }
-
-        // Add some empty polynomials to make Prover constructor happy; not used in the z_perm calculation
-        proving_key->polynomial_store.put("id_1_lagrange", polynomial(proving_key->circuit_size));
-        proving_key->polynomial_store.put("id_2_lagrange", polynomial(proving_key->circuit_size));
-        proving_key->polynomial_store.put("id_3_lagrange", polynomial(proving_key->circuit_size));
-        proving_key->polynomial_store.put("q_1_lagrange", polynomial(proving_key->circuit_size));
-        proving_key->polynomial_store.put("q_2_lagrange", polynomial(proving_key->circuit_size));
-        proving_key->polynomial_store.put("q_3_lagrange", polynomial(proving_key->circuit_size));
-        proving_key->polynomial_store.put("q_m_lagrange", polynomial(proving_key->circuit_size));
-        proving_key->polynomial_store.put("q_c_lagrange", polynomial(proving_key->circuit_size));
-        proving_key->polynomial_store.put("L_first_lagrange", polynomial(proving_key->circuit_size));
-        proving_key->polynomial_store.put("L_last_lagrange", polynomial(proving_key->circuit_size));
-
-        // Instantiate a Prover with wire polynomials and pointer to the proving_key just constructed
-        auto wires_copy = wires;
-        auto honk_prover = StandardProver(std::move(wires_copy), proving_key);
-
-        // Get random challenges
-        auto beta = Fscalar::random_element();
-        auto gamma = Fscalar::random_element();
-
-        // Method 1: Compute z_perm using 'compute_grand_product_polynomial' as the prover would in practice
-        polynomial prover_z_perm = honk_prover.compute_grand_product_polynomial(beta, gamma);
-
-        // Method 2: Compute z_perm locally using the simplest non-optimized syntax possible. The comment below,
-        // which describes the computation in 4 steps, is adapted from a similar comment in
-        // compute_grand_product_polynomial.
-        /*
-         * Assume program_width 3. Z_perm may be defined in terms of its values
-         * on X_i = 0,1,...,n-1 as Z_perm[0] = 0 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 three 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)
-         */
-
-        // Make scratch space for the numerator and denominator accumulators.
-        std::array<std::array<Fscalar, num_gates>, program_width> numererator_accum;
-        std::array<std::array<Fscalar, num_gates>, program_width> denominator_accum;
-
-        // Step (1)
-        for (size_t i = 0; i < proving_key->circuit_size; ++i) {
-            for (size_t k = 0; k < program_width; ++k) {
-                Fscalar idx = k * proving_key->circuit_size + i;                       // id_k[i]
-                numererator_accum[k][i] = wires[k][i] + (idx * beta) + gamma;          // w_k(i) + β.(k*n+i) + γ
-                denominator_accum[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 < proving_key->circuit_size - 1; ++i) {
-                numererator_accum[k][i + 1] *= numererator_accum[k][i];
-                denominator_accum[k][i + 1] *= denominator_accum[k][i];
-            }
-        }
-
-        // Step (3)
-        for (size_t i = 0; i < proving_key->circuit_size; ++i) {
-            for (size_t k = 1; k < program_width; ++k) {
-                numererator_accum[0][i] *= numererator_accum[k][i];
-                denominator_accum[0][i] *= denominator_accum[k][i];
-            }
-        }
-
-        // Step (4)
-        polynomial z_perm(proving_key->circuit_size);
-        z_perm[0] = Fscalar::zero(); // Z_0 = 1
-        // Note: in practice, we replace this expensive element-wise division with Montgomery batch inversion
-        for (size_t i = 0; i < proving_key->circuit_size - 1; ++i) {
-            z_perm[i + 1] = numererator_accum[0][i] / denominator_accum[0][i];
-        }
-
-        // Check consistency between locally computed z_perm and the one computed by the prover
-        EXPECT_EQ(z_perm, prover_z_perm);
-    };
-};
-
-typedef testing::Types<barretenberg::fr> FieldTypes;
-TYPED_TEST_SUITE(ProverTests, FieldTypes);
-
-TYPED_TEST(ProverTests, grand_product_construction)
-{
-    TestFixture::test_grand_product_construction();
-}
-
-} // namespace honk_prover_tests
diff --git a/cpp/src/barretenberg/honk/proof_system/prover_library.cpp b/cpp/src/barretenberg/honk/proof_system/prover_library.cpp
new file mode 100644
index 0000000000..fe102bcd2c
--- /dev/null
+++ b/cpp/src/barretenberg/honk/proof_system/prover_library.cpp
@@ -0,0 +1,360 @@
+#include "prover_library.hpp"
+#include "barretenberg/plonk/proof_system/types/prover_settings.hpp"
+#include <span>
+
+namespace honk::prover_library {
+
+using Fr = barretenberg::fr;
+using Polynomial = barretenberg::Polynomial<Fr>;
+
+/**
+ * @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 <size_t program_width>
+Polynomial compute_permutation_grand_product(std::shared_ptr<bonk::proving_key>& key,
+                                             std::vector<Polynomial>& wire_polynomials,
+                                             Fr beta,
+                                             Fr gamma)
+{
+    using barretenberg::polynomial_arithmetic::copy_polynomial;
+
+    // TODO(luke): instantiate z_perm here then make the 0th accum a span of it? avoid extra memory.
+
+    // Allocate accumulator polynomials that will serve as scratch space
+    std::array<Polynomial, program_width> numerator_accumulator;
+    std::array<Polynomial, program_width> denominator_accumulator;
+    for (size_t i = 0; i < program_width; ++i) {
+        numerator_accumulator[i] = Polynomial{ key->circuit_size };
+        denominator_accumulator[i] = Polynomial{ 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].
+    // Note: Since numerator_accumulator[0][i] corresponds to z_lookup[i+1], we only iterate up to i = (n - 2).
+    Fr* inversion_coefficients = &denominator_accumulator[1][0]; // arbitrary scratch space
+    Fr inversion_accumulator = Fr::one();
+    for (size_t i = 0; i < key->circuit_size - 1; ++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 - 2; 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);
+    // Initialize 0th coefficient to 0 to ensure z_perm is left-shiftable via division by X in gemini
+    z_perm[0] = 0;
+    copy_polynomial(numerator_accumulator[0].data(), &z_perm[1], key->circuit_size - 1, key->circuit_size - 1);
+
+    return z_perm;
+}
+
+/**
+ * @brief Compute the lookup grand product polynomial Z_lookup(X).
+ *
+ * @details The lookup grand product polynomial Z_lookup is of the form
+ *
+ *                   ∏(1 + β) ⋅ ∏(q_lookup*f_k + γ) ⋅ ∏(t_k + βt_{k+1} + γ(1 + β))
+ * Z_lookup(X_j) = -----------------------------------------------------------------
+ *                                   ∏(s_k + βs_{k+1} + γ(1 + β))
+ *
+ * where ∏ := ∏_{k<j}. This polynomial is constructed in evaluation form over the course
+ * of three steps:
+ *
+ * Step 1) Compute polynomials f, t and s and incorporate them into terms that are ultimately needed
+ * to construct the grand product polynomial Z_lookup(X):
+ * Note 1: In what follows, 't' is associated with table values (and is not to be confused with the
+ * quotient polynomial, also refered to as 't' elsewhere). Polynomial 's' is the sorted  concatenation
+ * of the witnesses and the table values.
+ * Note 2: Evaluation at Xω is indicated explicitly, e.g. 'p(Xω)'; evaluation at X is simply omitted, e.g. 'p'
+ *
+ * 1a.   Compute f, then set accumulators[0] = (q_lookup*f + γ), where
+ *
+ *         f = (w_1 + q_2*w_1(Xω)) + η(w_2 + q_m*w_2(Xω)) + η²(w_3 + q_c*w_3(Xω)) + η³q_index.
+ *      Note that q_2, q_m, and q_c are just the selectors from Standard Plonk that have been repurposed
+ *      in the context of the plookup gate to represent 'shift' values. For example, setting each of the
+ *      q_* in f to 2^8 facilitates operations on 32-bit values via four operations on 8-bit values. See
+ *      Ultra documentation for details.
+ *
+ * 1b.   Compute t, then set accumulators[1] = (t + βt(Xω) + γ(1 + β)), where t = t_1 + ηt_2 + η²t_3 + η³t_4
+ *
+ * 1c.   Set accumulators[2] = (1 + β)
+ *
+ * 1d.   Compute s, then set accumulators[3] = (s + βs(Xω) + γ(1 + β)), where s = s_1 + ηs_2 + η²s_3 + η³s_4
+ *
+ * Step 2) Compute the constituent product components of Z_lookup(X).
+ * Let ∏ := Prod_{k<j}. Let f_k, t_k and s_k now represent the k'th component of the polynomials f,t and s
+ * defined above. We compute the following four product polynomials needed to construct the grand product
+ * Z_lookup(X).
+ * 1.   accumulators[0][j] = ∏ (q_lookup*f_k + γ)
+ * 2.   accumulators[1][j] = ∏ (t_k + βt_{k+1} + γ(1 + β))
+ * 3.   accumulators[2][j] = ∏ (1 + β)
+ * 4.   accumulators[3][j] = ∏ (s_k + βs_{k+1} + γ(1 + β))
+ *
+ * Step 3) Combine the accumulator product elements to construct Z_lookup(X).
+ *
+ *                      ∏ (1 + β) ⋅ ∏ (q_lookup*f_k + γ) ⋅ ∏ (t_k + βt_{k+1} + γ(1 + β))
+ *  Z_lookup(g^j) = --------------------------------------------------------------------------
+ *                                      ∏ (s_k + βs_{k+1} + γ(1 + β))
+ *
+ * Note: Montgomery batch inversion is used to efficiently compute the coefficients of Z_lookup
+ * rather than peforming n individual inversions. I.e. we first compute the double product P_n:
+ *
+ * P_n := ∏_{j<n} ∏_{k<j} S_k, where S_k = (s_k + βs_{k+1} + γ(1 + β))
+ *
+ * and then compute the inverse on P_n. Then we work back to front to obtain terms of the form
+ * 1/∏_{k<i} S_i that appear in Z_lookup, using the fact that P_i/P_{i+1} = 1/∏_{k<i} S_i. (Note
+ * that once we have 1/P_n, we can compute 1/P_{n-1} as (1/P_n) * ∏_{k<n} S_i, and
+ * so on).
+ *
+ * @param key proving key
+ * @param wire_polynomials
+ * @param sorted_list_accumulator
+ * @param eta
+ * @param beta
+ * @param gamma
+ * @return Polynomial
+ */
+Polynomial compute_lookup_grand_product(std::shared_ptr<bonk::proving_key>& key,
+                                        std::vector<Polynomial>& wire_polynomials,
+                                        Polynomial& sorted_list_accumulator,
+                                        Fr eta,
+                                        Fr beta,
+                                        Fr gamma)
+{
+    const Fr eta_sqr = eta.sqr();
+    const Fr eta_cube = eta_sqr * eta;
+
+    const size_t circuit_size = key->circuit_size;
+
+    // Allocate 4 length n 'accumulator' polynomials. accumulators[0] will be used to construct
+    // z_lookup in place.
+    // Note: The magic number 4 here comes from the structure of the grand product and is not related to the program
+    // width.
+    std::array<Polynomial, 4> accumulators;
+    for (size_t i = 0; i < 4; ++i) {
+        accumulators[i] = Polynomial{ key->circuit_size };
+    }
+
+    // Obtain column step size values that have been stored in repurposed selctors
+    std::span<const Fr> column_1_step_size = key->polynomial_store.get("q_2_lagrange");
+    std::span<const Fr> column_2_step_size = key->polynomial_store.get("q_m_lagrange");
+    std::span<const Fr> column_3_step_size = key->polynomial_store.get("q_c_lagrange");
+
+    // Utilize three wires; this is not tied to program width
+    std::array<std::span<const Fr>, 3> wires;
+    for (size_t i = 0; i < 3; ++i) {
+        wires[i] = wire_polynomials[i];
+    }
+
+    // Note: the number of table polys is related to program width but '4' is the only value supported
+    std::array<std::span<const Fr>, 4> tables{
+        key->polynomial_store.get("table_value_1_lagrange"),
+        key->polynomial_store.get("table_value_2_lagrange"),
+        key->polynomial_store.get("table_value_3_lagrange"),
+        key->polynomial_store.get("table_value_4_lagrange"),
+    };
+
+    std::span<const Fr> lookup_selector = key->polynomial_store.get("table_type_lagrange");
+    std::span<const Fr> lookup_index_selector = key->polynomial_store.get("q_3_lagrange");
+
+    const Fr beta_plus_one = beta + Fr(1);                      // (1 + β)
+    const Fr gamma_times_beta_plus_one = gamma * beta_plus_one; // γ(1 + β)
+
+    // Step (1)
+
+    Fr T0; // intermediate value for various calculations below
+    // Note: block_mask is used for efficient modulus, i.e. i % N := i & (N-1), for N = 2^k
+    const size_t block_mask = circuit_size - 1;
+    // Initialize 't(X)' to be used in an expression of the form t(X) + β*t(Xω)
+    Fr next_table = tables[0][0] + tables[1][0] * eta + tables[2][0] * eta_sqr + tables[3][0] * eta_cube;
+
+    for (size_t i = 0; i < circuit_size; ++i) {
+
+        // Compute i'th element of f via Horner (see definition of f above)
+        T0 = lookup_index_selector[i];
+        T0 *= eta;
+        T0 += wires[2][(i + 1) & block_mask] * column_3_step_size[i];
+        T0 += wires[2][i];
+        T0 *= eta;
+        T0 += wires[1][(i + 1) & block_mask] * column_2_step_size[i];
+        T0 += wires[1][i];
+        T0 *= eta;
+        T0 += wires[0][(i + 1) & block_mask] * column_1_step_size[i];
+        T0 += wires[0][i];
+        T0 *= lookup_selector[i];
+
+        // Set i'th element of polynomial q_lookup*f + γ
+        accumulators[0][i] = T0;
+        accumulators[0][i] += gamma;
+
+        // Compute i'th element of t via Horner
+        T0 = tables[3][(i + 1) & block_mask];
+        T0 *= eta;
+        T0 += tables[2][(i + 1) & block_mask];
+        T0 *= eta;
+        T0 += tables[1][(i + 1) & block_mask];
+        T0 *= eta;
+        T0 += tables[0][(i + 1) & block_mask];
+
+        // Set i'th element of polynomial (t + βt(Xω) + γ(1 + β))
+        accumulators[1][i] = T0 * beta + next_table;
+        next_table = T0;
+        accumulators[1][i] += gamma_times_beta_plus_one;
+
+        // Set value of this accumulator to (1 + β)
+        accumulators[2][i] = beta_plus_one;
+
+        // Set i'th element of polynomial (s + βs(Xω) + γ(1 + β))
+        accumulators[3][i] = sorted_list_accumulator[(i + 1) & block_mask];
+        accumulators[3][i] *= beta;
+        accumulators[3][i] += sorted_list_accumulator[i];
+        accumulators[3][i] += gamma_times_beta_plus_one;
+    }
+
+    // Step (2)
+
+    // Note: This is a small multithreading bottleneck, as we have only 4 parallelizable processes.
+    for (auto& accum : accumulators) {
+        for (size_t i = 0; i < circuit_size - 1; ++i) {
+            accum[i + 1] *= accum[i];
+        }
+    }
+
+    // Step (3)
+
+    // Compute <Z_lookup numerator> * ∏_{j<i}∏_{k<j}S_k
+    Fr inversion_accumulator = Fr::one();
+    for (size_t i = 0; i < circuit_size - 1; ++i) {
+        accumulators[0][i] *= accumulators[2][i];
+        accumulators[0][i] *= accumulators[1][i];
+        accumulators[0][i] *= inversion_accumulator;
+        inversion_accumulator *= accumulators[3][i];
+    }
+    inversion_accumulator = inversion_accumulator.invert(); // invert
+
+    // Compute [Z_lookup numerator] * ∏_{j<i}∏_{k<j}S_k / ∏_{j<i+1}∏_{k<j}S_k = <Z_lookup numerator> /
+    // ∏_{k<i}S_k
+    for (size_t i = circuit_size - 2; i != std::numeric_limits<size_t>::max(); --i) {
+        // N.B. accumulators[0][i] = z_lookup[i + 1]
+        accumulators[0][i] *= inversion_accumulator;
+        inversion_accumulator *= accumulators[3][i];
+    }
+
+    Polynomial z_lookup(key->circuit_size);
+    // Initialize 0th coefficient to 0 to ensure z_perm is left-shiftable via division by X in gemini
+    z_lookup[0] = Fr::zero();
+    barretenberg::polynomial_arithmetic::copy_polynomial(
+        accumulators[0].data(), &z_lookup[1], key->circuit_size - 1, key->circuit_size - 1);
+
+    return z_lookup;
+}
+
+/**
+ * @brief Construct sorted list accumulator polynomial 's'.
+ *
+ * @details Compute s = s_1 + η*s_2 + η²*s_3 + η³*s_4 (via Horner) where s_i are the
+ * sorted concatenated witness/table polynomials
+ *
+ * @param key proving key
+ * @param sorted_list_polynomials sorted concatenated witness/table polynomials
+ * @param eta random challenge
+ * @return Polynomial
+ */
+Polynomial compute_sorted_list_accumulator(std::shared_ptr<bonk::proving_key>& key,
+                                           std::vector<Polynomial>& sorted_list_polynomials,
+                                           Fr eta)
+{
+    const size_t circuit_size = key->circuit_size;
+
+    barretenberg::polynomial sorted_list_accumulator(sorted_list_polynomials[0]);
+    std::span<const Fr> s_2 = sorted_list_polynomials[1];
+    std::span<const Fr> s_3 = sorted_list_polynomials[2];
+    std::span<const Fr> s_4 = sorted_list_polynomials[3];
+
+    // Construct s via Horner, i.e. s = s_1 + η(s_2 + η(s_3 + η*s_4))
+    for (size_t i = 0; i < circuit_size; ++i) {
+        Fr T0 = s_4[i];
+        T0 *= eta;
+        T0 += s_3[i];
+        T0 *= eta;
+        T0 += s_2[i];
+        T0 *= eta;
+        sorted_list_accumulator[i] += T0;
+    }
+
+    return sorted_list_accumulator;
+}
+
+template Polynomial compute_permutation_grand_product<plonk::standard_settings::program_width>(
+    std::shared_ptr<bonk::proving_key>&, std::vector<Polynomial>&, Fr, Fr);
+template Polynomial compute_permutation_grand_product<plonk::ultra_settings::program_width>(
+    std::shared_ptr<bonk::proving_key>&, std::vector<Polynomial>&, Fr, Fr);
+
+} // namespace honk::prover_library
diff --git a/cpp/src/barretenberg/honk/proof_system/prover_library.hpp b/cpp/src/barretenberg/honk/proof_system/prover_library.hpp
new file mode 100644
index 0000000000..0fc6d1662f
--- /dev/null
+++ b/cpp/src/barretenberg/honk/proof_system/prover_library.hpp
@@ -0,0 +1,30 @@
+#pragma once
+#include "barretenberg/ecc/curves/bn254/fr.hpp"
+#include "barretenberg/polynomials/polynomial.hpp"
+#include "barretenberg/proof_system/proving_key/proving_key.hpp"
+#include "barretenberg/plonk/proof_system/types/proof.hpp"
+#include "barretenberg/plonk/proof_system/types/program_settings.hpp"
+
+namespace honk::prover_library {
+
+using Fr = barretenberg::fr;
+using Polynomial = barretenberg::Polynomial<Fr>;
+
+template <size_t program_width>
+Polynomial compute_permutation_grand_product(std::shared_ptr<bonk::proving_key>& key,
+                                             std::vector<Polynomial>& wire_polynomials,
+                                             Fr beta,
+                                             Fr gamma);
+
+Polynomial compute_lookup_grand_product(std::shared_ptr<bonk::proving_key>& key,
+                                        std::vector<Polynomial>& wire_polynomials,
+                                        Polynomial& sorted_list_accumulator,
+                                        Fr eta,
+                                        Fr beta,
+                                        Fr gamma);
+
+Polynomial compute_sorted_list_accumulator(std::shared_ptr<bonk::proving_key>& key,
+                                           std::vector<Polynomial>& sorted_list_polynomials,
+                                           Fr eta);
+
+} // namespace honk::prover_library
diff --git a/cpp/src/barretenberg/honk/proof_system/prover_library.test.cpp b/cpp/src/barretenberg/honk/proof_system/prover_library.test.cpp
new file mode 100644
index 0000000000..205d249581
--- /dev/null
+++ b/cpp/src/barretenberg/honk/proof_system/prover_library.test.cpp
@@ -0,0 +1,340 @@
+#include "barretenberg/plonk/proof_system/types/prover_settings.hpp"
+#include "prover.hpp"
+#include "prover_library.hpp"
+#include "barretenberg/polynomials/polynomial.hpp"
+
+#include "barretenberg/srs/reference_string/file_reference_string.hpp"
+#include <array>
+#include <string>
+#include <vector>
+#include <cstddef>
+#include <gtest/gtest.h>
+
+using namespace honk;
+namespace prover_library_tests {
+
+// field is named Fscalar here because of clash with the Fr
+template <class FF> class ProverLibraryTests : public testing::Test {
+
+    using Polynomial = barretenberg::Polynomial<FF>;
+
+  public:
+    /**
+     * @brief Get a random polynomial
+     *
+     * @param size
+     * @return Polynomial
+     */
+    static constexpr Polynomial get_random_polynomial(size_t size)
+    {
+        Polynomial random_polynomial{ size };
+        for (auto& coeff : random_polynomial) {
+            coeff = FF::random_element();
+        }
+        return random_polynomial;
+    }
+
+    /**
+     * @brief Check consistency of the computation of the permutation grand product polynomial z_permutation.
+     * @details This test compares a simple, unoptimized, easily readable calculation of the grand product z_permutation
+     * to the optimized implementation used by the prover. It's purpose is to provide confidence that some optimization
+     * introduced into the calculation has not changed the result.
+     * @note This test does confirm the correctness of z_permutation, only that the two implementations yield an
+     * identical result.
+     */
+    template <size_t program_width> static void test_permutation_grand_product_construction()
+    {
+        // Define some mock inputs for proving key constructor
+        static const size_t num_gates = 8;
+        static const size_t num_public_inputs = 0;
+        auto reference_string = std::make_shared<bonk::FileReferenceString>(num_gates + 1, "../srs_db/ignition");
+
+        // Instatiate a proving_key and make a pointer to it. This will be used to instantiate a Prover.
+        auto proving_key = std::make_shared<bonk::proving_key>(
+            num_gates, num_public_inputs, reference_string, plonk::ComposerType::STANDARD_HONK);
+
+        // static const size_t program_width = StandardProver::settings_::program_width;
+
+        // Construct mock wire and permutation polynomials.
+        // Note: for the purpose of checking the consistency between two methods of computing z_perm, these polynomials
+        // can simply be random. We're not interested in the particular properties of the result.
+        std::vector<Polynomial> wires;
+        std::vector<Polynomial> sigmas;
+        for (size_t i = 0; i < program_width; ++i) {
+            wires.emplace_back(get_random_polynomial(num_gates));
+            sigmas.emplace_back(get_random_polynomial(num_gates));
+
+            // Add sigma polys to proving_key; to be used by the prover in constructing it's own z_perm
+            std::string sigma_id = "sigma_" + std::to_string(i + 1) + "_lagrange";
+            proving_key->polynomial_store.put(sigma_id, Polynomial{ sigmas[i] });
+        }
+
+        // Get random challenges
+        auto beta = FF::random_element();
+        auto gamma = FF::random_element();
+
+        // Method 1: Compute z_perm using 'compute_grand_product_polynomial' as the prover would in practice
+        Polynomial z_permutation =
+            prover_library::compute_permutation_grand_product<program_width>(proving_key, wires, beta, gamma);
+
+        // Method 2: Compute z_perm locally using the simplest non-optimized syntax possible. The comment below,
+        // which describes the computation in 4 steps, is adapted from a similar comment in
+        // compute_grand_product_polynomial.
+        /*
+         * Assume program_width 3. Z_perm may be defined in terms of its values
+         * on X_i = 0,1,...,n-1 as Z_perm[0] = 0 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 three 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)
+         */
+
+        // Make scratch space for the numerator and denominator accumulators.
+        std::array<std::array<FF, num_gates>, program_width> numererator_accum;
+        std::array<std::array<FF, num_gates>, program_width> denominator_accum;
+
+        // Step (1)
+        for (size_t i = 0; i < proving_key->circuit_size; ++i) {
+            for (size_t k = 0; k < program_width; ++k) {
+                FF idx = k * proving_key->circuit_size + i;                            // id_k[i]
+                numererator_accum[k][i] = wires[k][i] + (idx * beta) + gamma;          // w_k(i) + β.(k*n+i) + γ
+                denominator_accum[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 < proving_key->circuit_size - 1; ++i) {
+                numererator_accum[k][i + 1] *= numererator_accum[k][i];
+                denominator_accum[k][i + 1] *= denominator_accum[k][i];
+            }
+        }
+
+        // Step (3)
+        for (size_t i = 0; i < proving_key->circuit_size; ++i) {
+            for (size_t k = 1; k < program_width; ++k) {
+                numererator_accum[0][i] *= numererator_accum[k][i];
+                denominator_accum[0][i] *= denominator_accum[k][i];
+            }
+        }
+
+        // Step (4)
+        Polynomial z_permutation_expected(proving_key->circuit_size);
+        z_permutation_expected[0] = FF::zero(); // Z_0 = 1
+        // Note: in practice, we replace this expensive element-wise division with Montgomery batch inversion
+        for (size_t i = 0; i < proving_key->circuit_size - 1; ++i) {
+            z_permutation_expected[i + 1] = numererator_accum[0][i] / denominator_accum[0][i];
+        }
+
+        // Check consistency between locally computed z_perm and the one computed by the prover library
+        EXPECT_EQ(z_permutation, z_permutation_expected);
+    };
+
+    /**
+     * @brief Check consistency of the computation of the lookup grand product polynomial z_lookup.
+     * @details This test compares a simple, unoptimized, easily readable calculation of the grand product z_lookup
+     * to the optimized implementation used by the prover. It's purpose is to provide confidence that some optimization
+     * introduced into the calculation has not changed the result.
+     * @note This test does confirm the correctness of z_lookup, only that the two implementations yield an
+     * identical result.
+     */
+    static void test_lookup_grand_product_construction()
+    {
+        // Define some mock inputs for proving key constructor
+        static const size_t circuit_size = 8;
+        static const size_t num_public_inputs = 0;
+        auto reference_string = std::make_shared<bonk::FileReferenceString>(circuit_size + 1, "../srs_db/ignition");
+
+        // Instatiate a proving_key and make a pointer to it. This will be used to instantiate a Prover.
+        auto proving_key = std::make_shared<bonk::proving_key>(
+            circuit_size, num_public_inputs, reference_string, plonk::ComposerType::STANDARD_HONK);
+
+        // Construct mock wire and permutation polynomials.
+        // Note: for the purpose of checking the consistency between two methods of computing z_perm, these polynomials
+        // can simply be random. We're not interested in the particular properties of the result.
+        std::vector<Polynomial> wires;
+        for (size_t i = 0; i < 3; ++i) {
+            wires.emplace_back(get_random_polynomial(circuit_size));
+        }
+        std::vector<Polynomial> tables;
+        for (size_t i = 0; i < 4; ++i) {
+            tables.emplace_back(get_random_polynomial(circuit_size));
+            std::string label = "table_value_" + std::to_string(i + 1) + "_lagrange";
+            proving_key->polynomial_store.put(label, Polynomial{ tables[i] });
+        }
+
+        auto s_lagrange = get_random_polynomial(circuit_size);
+        auto column_1_step_size = get_random_polynomial(circuit_size);
+        auto column_2_step_size = get_random_polynomial(circuit_size);
+        auto column_3_step_size = get_random_polynomial(circuit_size);
+        auto lookup_index_selector = get_random_polynomial(circuit_size);
+        auto lookup_selector = get_random_polynomial(circuit_size);
+
+        proving_key->polynomial_store.put("s_lagrange", Polynomial{ s_lagrange });
+        proving_key->polynomial_store.put("q_2_lagrange", Polynomial{ column_1_step_size });
+        proving_key->polynomial_store.put("q_m_lagrange", Polynomial{ column_2_step_size });
+        proving_key->polynomial_store.put("q_c_lagrange", Polynomial{ column_3_step_size });
+        proving_key->polynomial_store.put("q_3_lagrange", Polynomial{ lookup_index_selector });
+        proving_key->polynomial_store.put("table_type_lagrange", Polynomial{ lookup_selector });
+
+        // Get random challenges
+        auto beta = FF::random_element();
+        auto gamma = FF::random_element();
+        auto eta = FF::random_element();
+
+        // Method 1: Compute z_lookup using the prover library method
+        Polynomial z_lookup =
+            prover_library::compute_lookup_grand_product(proving_key, wires, s_lagrange, eta, beta, gamma);
+
+        // Method 2: Compute the lookup grand product polynomial Z_lookup:
+        //
+        //                   ∏(1 + β) ⋅ ∏(q_lookup*f_k + γ) ⋅ ∏(t_k + βt_{k+1} + γ(1 + β))
+        // Z_lookup(X_j) = -----------------------------------------------------------------
+        //                                   ∏(s_k + βs_{k+1} + γ(1 + β))
+        //
+        // in a way that is simple to read (but inefficient). See prover library method for more details.
+        const Fr eta_sqr = eta.sqr();
+        const Fr eta_cube = eta_sqr * eta;
+
+        std::array<Polynomial, 4> accumulators;
+        for (size_t i = 0; i < 4; ++i) {
+            accumulators[i] = Polynomial{ circuit_size };
+        }
+
+        // Step (1)
+
+        // Note: block_mask is used for efficient modulus, i.e. i % N := i & (N-1), for N = 2^k
+        const size_t block_mask = circuit_size - 1;
+        // Initialize 't(X)' to be used in an expression of the form t(X) + β*t(Xω)
+        Fr table_i = tables[0][0] + tables[1][0] * eta + tables[2][0] * eta_sqr + tables[3][0] * eta_cube;
+        for (size_t i = 0; i < circuit_size; ++i) {
+            size_t shift_idx = (i + 1) & block_mask;
+
+            // f = (w_1 + q_2*w_1(Xω)) + η(w_2 + q_m*w_2(Xω)) + η²(w_3 + q_c*w_3(Xω)) + η³q_index.
+            Fr f_i = (wires[0][i] + wires[0][shift_idx] * column_1_step_size[i]) +
+                     (wires[1][i] + wires[1][shift_idx] * column_2_step_size[i]) * eta +
+                     (wires[2][i] + wires[2][shift_idx] * column_3_step_size[i]) * eta_sqr +
+                     eta_cube * lookup_index_selector[i];
+
+            // q_lookup * f + γ
+            accumulators[0][i] = lookup_selector[i] * f_i + gamma;
+
+            // t = t_1 + ηt_2 + η²t_3 + η³t_4
+            Fr table_i_plus_1 = tables[0][shift_idx] + eta * tables[1][shift_idx] + eta_sqr * tables[2][shift_idx] +
+                                eta_cube * tables[3][shift_idx];
+
+            // t + βt(Xω) + γ(1 + β)
+            accumulators[1][i] = table_i + table_i_plus_1 * beta + gamma * (Fr::one() + beta);
+
+            // (1 + β)
+            accumulators[2][i] = Fr::one() + beta;
+
+            // s + βs(Xω) + γ(1 + β)
+            accumulators[3][i] = s_lagrange[i] + beta * s_lagrange[shift_idx] + gamma * (Fr::one() + beta);
+
+            // Set t(X_i) for next iteration
+            table_i = table_i_plus_1;
+        }
+
+        // Step (2)
+        for (auto& accum : accumulators) {
+            for (size_t i = 0; i < circuit_size - 1; ++i) {
+                accum[i + 1] *= accum[i];
+            }
+        }
+
+        // Step (3)
+        Polynomial z_lookup_expected(circuit_size);
+        z_lookup_expected[0] = FF::zero(); // Z_lookup_0 = 0
+
+        // Compute the numerator in accumulators[0]; The denominator is in accumulators[3]
+        for (size_t i = 0; i < circuit_size - 1; ++i) {
+            accumulators[0][i] *= accumulators[1][i] * accumulators[2][i];
+        }
+        // Compute Z_lookup_i, i = [1, n-1]
+        for (size_t i = 0; i < circuit_size - 1; ++i) {
+            z_lookup_expected[i + 1] = accumulators[0][i] / accumulators[3][i];
+        }
+
+        EXPECT_EQ(z_lookup, z_lookup_expected);
+    };
+
+    /**
+     * @brief Check consistency of the computation of the sorted list accumulator
+     * @details This test compares a simple, unoptimized, easily readable calculation of the sorted list accumulator
+     * to the optimized implementation used by the prover. It's purpose is to provide confidence that some optimization
+     * introduced into the calculation has not changed the result.
+     * @note This test does confirm the correctness of the sorted list accumulator, only that the two implementations
+     * yield an identical result.
+     */
+    static void test_sorted_list_accumulator_construction()
+    {
+        // Construct a proving_key
+        static const size_t circuit_size = 8;
+        static const size_t num_public_inputs = 0;
+        auto reference_string = std::make_shared<bonk::FileReferenceString>(circuit_size + 1, "../srs_db/ignition");
+        auto proving_key = std::make_shared<bonk::proving_key>(
+            circuit_size, num_public_inputs, reference_string, plonk::ComposerType::STANDARD_HONK);
+
+        // Get random challenge eta
+        auto eta = FF::random_element();
+
+        // Construct mock sorted list polynomials.
+        std::vector<Polynomial> sorted_list_polynomials;
+        for (size_t i = 0; i < 4; ++i) {
+            sorted_list_polynomials.emplace_back(get_random_polynomial(circuit_size));
+        }
+
+        // Method 1: computed sorted list accumulator polynomial using prover library method
+        Polynomial sorted_list_accumulator =
+            prover_library::compute_sorted_list_accumulator(proving_key, sorted_list_polynomials, eta);
+
+        // Method 2: Compute local sorted list accumulator simply and inefficiently
+        const Fr eta_sqr = eta.sqr();
+        const Fr eta_cube = eta_sqr * eta;
+
+        // Compute s = s_1 + η*s_2 + η²*s_3 + η³*s_4
+        Polynomial sorted_list_accumulator_expected{ sorted_list_polynomials[0] };
+        for (size_t i = 0; i < circuit_size; ++i) {
+            sorted_list_accumulator_expected[i] += sorted_list_polynomials[1][i] * eta +
+                                                   sorted_list_polynomials[2][i] * eta_sqr +
+                                                   sorted_list_polynomials[3][i] * eta_cube;
+        }
+
+        EXPECT_EQ(sorted_list_accumulator, sorted_list_accumulator_expected);
+    };
+};
+
+typedef testing::Types<barretenberg::fr> FieldTypes;
+TYPED_TEST_SUITE(ProverLibraryTests, FieldTypes);
+
+TYPED_TEST(ProverLibraryTests, PermutationGrandProduct)
+{
+    TestFixture::template test_permutation_grand_product_construction<plonk::standard_settings::program_width>();
+    TestFixture::template test_permutation_grand_product_construction<plonk::ultra_settings::program_width>();
+}
+
+TYPED_TEST(ProverLibraryTests, LookupGrandProduct)
+{
+    TestFixture::test_lookup_grand_product_construction();
+}
+
+TYPED_TEST(ProverLibraryTests, SortedListAccumulator)
+{
+    TestFixture::test_sorted_list_accumulator_construction();
+}
+
+} // namespace prover_library_tests