barycentric.hpp

Go to the documentation of this file.

// === AUDIT STATUS ===
// internal:    { status: Planned, auditors: [], commit: }
// external_1:  { status: not started, auditors: [], commit: }
// external_2:  { status: not started, auditors: [], commit: }
// =====================
 
#pragma once
#include "barretenberg/common/assert.hpp"
#include "barretenberg/ecc/fields/field.hpp"
#include <array>
 
/* Future improvements (see https://github.com/AztecProtocol/barretenberg/issues/10): The code works for its intended
 * use but could be improved: Precomputing for all possible size pairs is
 * probably feasible and might be a better solution than instantiating many instances separately. Then perhaps we could
 * infer input type to `extend`.
 */
namespace bb {
 
template <typename T> struct is_field_type {
    static constexpr bool value = false;
};
 
template <typename Params> struct is_field_type<bb::field<Params>> {
    static constexpr bool value = true;
};
 
template <typename T> inline constexpr bool is_field_type_v = is_field_type<T>::value;
 
template <class Fr, size_t domain_end, size_t num_evals> class BarycentricDataFunctions {
  public:
    static constexpr size_t domain_size = domain_end;
    static constexpr size_t big_domain_size = std::max(domain_size, num_evals);
 
    static_assert(domain_size > 0, "BarycentricData requires domain_size > 0");
    static_assert(num_evals > 0, "BarycentricData requires num_evals > 0");
    static_assert(num_evals >= domain_end || (!is_field_type_v<Fr> && num_evals == 1),
                  "Expected num_evals >= domain_size (or stdlib num_evals==1 special case)");
 
    // build big_domain, currently the set of x_i in {0, ..., big_domain_size - 1 }
    static constexpr std::array<Fr, big_domain_size> construct_big_domain()
    {
        std::array<Fr, big_domain_size> result;
        for (size_t i = 0; i < big_domain_size; ++i) {
            result[i] = static_cast<Fr>(i);
        }
        return result;
    }
 
    // build set of lagrange_denominators d_i = \prod_{j!=i} x_i - x_j
    static constexpr std::array<Fr, domain_size> construct_lagrange_denominators(const auto& big_domain)
    {
        std::array<Fr, domain_size> result;
        for (size_t i = 0; i != domain_size; ++i) {
            result[i] = 1;
            for (size_t j = 0; j != domain_size; ++j) {
                if (j != i) {
                    result[i] *= big_domain[i] - big_domain[j];
                }
            }
        }
        return result;
    }
 
    // Non-constexpr helper so we can use BB_ASSERT inside the constexpr batch_invert.
    static void assert_constant(const Fr& val)
    {
        if constexpr (!is_field_type_v<Fr>) {
            BB_ASSERT(val.is_constant(), "barycentric coeffs must be constants, not witnesses");
        }
    }
 
    static constexpr std::array<Fr, domain_size * num_evals> batch_invert(
        const std::array<Fr, domain_size * num_evals>& coeffs)
    {
        constexpr size_t n = domain_size * num_evals;
        std::array<Fr, n> temporaries{};
        std::array<bool, n> skipped{};
        Fr accumulator = 1;
        for (size_t i = 0; i < n; ++i) {
            temporaries[i] = accumulator;
            // For native field types, == 0 is a plain constexpr comparison. For stdlib field_t,
            // operator== would create circuit constraints and return bool_t, so we use get_value()
            // to extract the underlying native value for a plain comparison instead. Note that coeffs[] are derived
            // from constants (big_domain[i] = Fr(i)) and products/differences thereof, so they are constants. This
            // makes the get_value() based zero check safe.
            bool is_zero = false;
            if constexpr (is_field_type_v<Fr>) {
                is_zero = (coeffs[i] == 0);
            } else {
                assert_constant(coeffs[i]);
                is_zero = (coeffs[i].get_value() == 0);
            }
            if (is_zero) {
                skipped[i] = true;
            } else {
                skipped[i] = false;
                accumulator *= coeffs[i];
            }
        }
        accumulator = Fr(1) / accumulator;
        std::array<Fr, n> result{};
        Fr T0;
        for (size_t i = n - 1; i < n; --i) {
            if (!skipped[i]) {
                T0 = accumulator * temporaries[i];
                accumulator *= coeffs[i];
                result[i] = T0;
            }
        }
        return result;
    }
 
    // for each x_k in the big domain, build set of domain size-many denominator inverses
    // 1/(d_i*(x_k - x_j)). will multiply against each of these (rather than to divide by something)
    // for each barycentric evaluation
    // special case for stdlib path: if num_evals == 1, we output the barycentric weights result[j] = 1 / d_j.
    // This case does not arise on the native (compile-time) path.
    static constexpr std::array<Fr, domain_size * num_evals> construct_denominator_inverses(
        const auto& big_domain, const auto& lagrange_denominators)
    {
        std::array<Fr, domain_size * num_evals> result{}; // default init to 0 since below does not init all elements
        if constexpr (!is_field_type_v<Fr> && num_evals == 1) {
            result = lagrange_denominators;
        } else {
            // Used in Univariate's `extend_to` method to extend univariates given by > 4 evaluations ( deg>3 ) to a
            // bigger evaluation domain.
            for (size_t k = domain_size; k < num_evals; ++k) {
                for (size_t j = 0; j < domain_size; ++j) {
                    Fr inv = lagrange_denominators[j];
                    inv *= (big_domain[k] - big_domain[j]);
                    result[(k * domain_size) + j] = inv;
                }
            }
        }
        return batch_invert(result);
    }
 
    // get full numerator values
    // full numerator is M(x) = \prod_{i} (x-x_i)
    // these will be zero for i < domain_size, but that's ok because
    // at such entries we will already have the evaluations of the polynomial
    static constexpr std::array<Fr, num_evals> construct_full_numerator_values(const auto& big_domain)
    {
        std::array<Fr, num_evals> result;
        for (size_t i = 0; i != num_evals; ++i) {
            result[i] = 1;
            Fr v_i = i;
            for (size_t j = 0; j != domain_size; ++j) {
                result[i] *= v_i - big_domain[j];
            }
        }
        return result;
    }
};
 
template <class Fr, size_t domain_end, size_t num_evals>
class BarycentricDataCompileTime : public BarycentricDataFunctions<Fr, domain_end, num_evals> {
    using Base = BarycentricDataFunctions<Fr, domain_end, num_evals>;
 
  public:
    using Base::big_domain_size;
    using Base::domain_size;
 
    static constexpr auto big_domain = Base::construct_big_domain();
    static constexpr auto lagrange_denominators = Base::construct_lagrange_denominators(big_domain);
    static constexpr auto precomputed_denominator_inverses =
        Base::construct_denominator_inverses(big_domain, lagrange_denominators);
    static constexpr auto full_numerator_values = Base::construct_full_numerator_values(big_domain);
};
 
template <class Fr, size_t domain_end, size_t num_evals>
class BarycentricDataRunTime : public BarycentricDataFunctions<Fr, domain_end, num_evals> {
    using Base = BarycentricDataFunctions<Fr, domain_end, num_evals>;
 
  public:
    using Base::big_domain_size;
    using Base::domain_size;
 
    inline static const auto big_domain = Base::construct_big_domain();
    inline static const auto lagrange_denominators = Base::construct_lagrange_denominators(big_domain);
    inline static const auto precomputed_denominator_inverses =
        Base::construct_denominator_inverses(big_domain, lagrange_denominators);
    inline static const auto full_numerator_values = Base::construct_full_numerator_values(big_domain);
};
 
template <class Fr, size_t domain_end, size_t num_evals>
using BarycentricData = std::conditional_t<is_field_type_v<Fr>,
                                           BarycentricDataCompileTime<Fr, domain_end, num_evals>,
                                           BarycentricDataRunTime<Fr, domain_end, num_evals>>;
 
} // namespace bb