grand_product_library.hpp

Go to the documentation of this file.

// === AUDIT STATUS ===
// internal:    { status: Complete, auditors: [Raju], commit: 21a7e3670e6 }
// external_1:  { status: not started, auditors: [], commit: }
// external_2:  { status: not started, auditors: [], commit: }
// =====================
 
#pragma once
#include "barretenberg/common/assert.hpp"
#include "barretenberg/common/constexpr_utils.hpp"
#include "barretenberg/common/debug_log.hpp"
#include "barretenberg/common/thread.hpp"
#include "barretenberg/common/zip_view.hpp"
#include "barretenberg/constants.hpp"
#include "barretenberg/flavor/flavor.hpp"
 
#include "barretenberg/relations/relation_parameters.hpp"
#include "barretenberg/trace_to_polynomials/trace_to_polynomials.hpp"
 
namespace bb {
 
template <typename Flavor, typename GrandProdRelation>
void compute_grand_product(typename Flavor::ProverPolynomials& full_polynomials,
                           bb::RelationParameters<typename Flavor::FF>& relation_parameters,
                           size_t size_override = 0,
                           uint32_t* duplicate_count_out = nullptr)
{
    BB_BENCH_NAME("compute_grand_product");
 
    using FF = typename Flavor::FF;
    using Polynomial = typename Flavor::Polynomial;
    using Accumulator = std::tuple_element_t<0, typename GrandProdRelation::SumcheckArrayOfValuesOverSubrelations>;
 
    // Set the domain over which the grand product must be computed. This may be less than the dyadic circuit size, e.g
    // the permutation grand product does not need to be computed beyond the index of the last active wire
    size_t domain_size = size_override == 0 ? full_polynomials.get_polynomial_size() : size_override;
 
    // Grand product starts after the disabled/reserved head rows (where lagrange_first lives).
    constexpr size_t gp_start = Flavor::TRACE_OFFSET;
 
    const size_t active_size = domain_size - gp_start;
 
    // The size of the iteration domain is one less than the active domain since the final value of the
    // grand product is constructed only in the relation and not explicitly in the polynomial
    const MultithreadData thread_data = calculate_thread_data(active_size - 1);
 
    // Allocate numerator/denominator polynomials that will serve as scratch space
    // TODO: we can re-use the permutation polynomial as the numerator polynomial (reduces readability)
    Polynomial numerator{ active_size };
    Polynomial denominator{ active_size };
 
    // Optionally count rows whose per-row ratio is 1 (numerator == denominator): there z_perm is
    // unchanged across the row, i.e. z_perm[i] == z_perm[i-1]. These are exactly the adjacent
    // duplicate coefficients the MSM dedup pre-pass strips. Counting them here (the per-row values
    // are already in hand) avoids a second full pass over z_perm in the caller.
    const bool count_duplicates = duplicate_count_out != nullptr;
    std::vector<uint32_t> thread_duplicate_counts(count_duplicates ? thread_data.num_threads : 0, 0);
 
    // Step (1)
    // Populate `numerator` and `denominator` with the algebra described by Relation
    parallel_for(thread_data.num_threads, [&](size_t thread_idx) {
        BB_BENCH_TRACY_NAME("GrandProduct::step1_numerator_denominator");
        const size_t start = thread_data.start[thread_idx];
        const size_t end = thread_data.end[thread_idx];
        typename Flavor::AllValues row;
        uint32_t local_duplicates = 0;
        for (size_t i = start; i < end; ++i) {
            const size_t poly_idx = i + gp_start;
            // TODO: consider avoiding get_row if possible.
            if constexpr (IsUltraOrMegaHonk<Flavor>) {
                row = full_polynomials.get_row_for_permutation_arg(poly_idx);
            } else {
                row = full_polynomials.get_row(poly_idx);
            }
            numerator.at(i) =
                GrandProdRelation::template compute_grand_product_numerator<Accumulator>(row, relation_parameters);
            denominator.at(i) =
                GrandProdRelation::template compute_grand_product_denominator<Accumulator>(row, relation_parameters);
            if (count_duplicates && numerator[i] == denominator[i]) {
                ++local_duplicates;
            }
        }
        if (count_duplicates) {
            thread_duplicate_counts[thread_idx] = local_duplicates;
        }
    });
 
    if (count_duplicates) {
        uint32_t total_duplicates = 0;
        for (const uint32_t c : thread_duplicate_counts) {
            total_duplicates += c;
        }
        *duplicate_count_out = total_duplicates;
    }
 
    DEBUG_LOG_ALL(numerator.coeffs());
    DEBUG_LOG_ALL(denominator.coeffs());
 
    // Step (2)
    // Compute the accumulating product of the numerator and denominator terms.
    // This step is split into three parts for efficient multithreading:
    // (i) compute ∏ A(j), ∏ B(j) subproducts for each thread
    // (ii) compute scaling factor required to convert each subproduct into a single running product
    // (ii) combine subproducts into a single running product
    //
    // For example, consider 4 threads and a size-8 numerator { a0, a1, a2, a3, a4, a5, a6, a7 }
    // (i)   Each thread computes 1 element of N = {{ a0, a0a1 }, { a2, a2a3 }, { a4, a4a5 }, { a6, a6a7 }}
    // (ii)  Take partial products P = { 1, a0a1, a2a3, a4a5 }
    // (iii) Each thread j computes N[i][j]*P[j]=
    //      {{a0,a0a1},{a0a1a2,a0a1a2a3},{a0a1a2a3a4,a0a1a2a3a4a5},{a0a1a2a3a4a5a6,a0a1a2a3a4a5a6a7}}
    std::vector<FF> partial_numerators(thread_data.num_threads);
    std::vector<FF> partial_denominators(thread_data.num_threads);
 
    parallel_for(thread_data.num_threads, [&](size_t thread_idx) {
        BB_BENCH_TRACY_NAME("GrandProduct::step2a_subproducts");
        const size_t start = thread_data.start[thread_idx];
        const size_t end = thread_data.end[thread_idx];
        for (size_t i = start; i < end - 1; ++i) {
            numerator.at(i + 1) *= numerator[i];
            denominator.at(i + 1) *= denominator[i];
        }
        partial_numerators[thread_idx] = numerator[end - 1];
        partial_denominators[thread_idx] = denominator[end - 1];
    });
 
    DEBUG_LOG_ALL(partial_numerators);
    DEBUG_LOG_ALL(partial_denominators);
 
    parallel_for(thread_data.num_threads, [&](size_t thread_idx) {
        BB_BENCH_TRACY_NAME("GrandProduct::step2b_scale_and_invert");
        const size_t start = thread_data.start[thread_idx];
        const size_t end = thread_data.end[thread_idx];
        if (thread_idx > 0) {
            FF numerator_scaling = 1;
            FF denominator_scaling = 1;
 
            for (size_t j = 0; j < thread_idx; ++j) {
                numerator_scaling *= partial_numerators[j];
                denominator_scaling *= partial_denominators[j];
            }
            for (size_t i = start; i < end; ++i) {
                numerator.at(i) = numerator[i] * numerator_scaling;
                denominator.at(i) = denominator[i] * denominator_scaling;
            }
        }
 
        // Final step: invert denominator
        FF::batch_invert(std::span{ &denominator.data()[start], end - start });
    });
 
    DEBUG_LOG_ALL(numerator.coeffs());
    DEBUG_LOG_ALL(denominator.coeffs());
 
    // Step (3) Compute grand_product_polynomial[i] = numerator[i] / denominator[i]
    auto& grand_product_polynomial = GrandProdRelation::get_grand_product_polynomial(full_polynomials);
    // The grand_product_polynomial must be shiftable for the permutation argument
    BB_ASSERT(grand_product_polynomial.is_shiftable());
 
    // Initialize grand product: z_perm[gp_start] = 1 (the first active row after disabled region)
    // For non-ZK, gp_start = NUM_ZERO_ROWS = 1, which matches z_perm[1] = num[0] * inv_den[0] implicitly
    // (since z_perm[0] = 0, the relation at lagrange_first uses z_perm + 1).
    // For ZK with top masking, z_perm must be 0 at lagrange_first (row NUM_DISABLED_ROWS_IN_SUMCHECK)
    // and the product starts from gp_start.
 
    // Compute grand product values: z_perm[gp_start + i + 1] = numerator[i] / denominator[i]
    parallel_for(thread_data.num_threads, [&](size_t thread_idx) {
        BB_BENCH_TRACY_NAME("GrandProduct::step3_quotient");
        const size_t start = thread_data.start[thread_idx];
        const size_t end = thread_data.end[thread_idx];
        for (size_t i = start; i < end; ++i) {
            grand_product_polynomial.at(gp_start + i + 1) = numerator[i] * denominator[i];
        }
    });
 
    DEBUG_LOG_ALL(grand_product_polynomial.coeffs());
}
 
template <typename Flavor>
void compute_grand_products(typename Flavor::ProverPolynomials& full_polynomials,
                            bb::RelationParameters<typename Flavor::FF>& relation_parameters,
                            const size_t size_override = 0)
{
    using GrandProductRelations = typename Flavor::GrandProductRelations;
 
    constexpr size_t NUM_RELATIONS = std::tuple_size<GrandProductRelations>{};
    bb::constexpr_for<0, NUM_RELATIONS, 1>([&]<size_t i>() {
        using GrandProdRelation = typename std::tuple_element<i, GrandProductRelations>::type;
 
        compute_grand_product<Flavor, GrandProdRelation>(full_polynomials, relation_parameters, size_override);
    });
}
 
} // namespace bb