univariate.hpp

Go to the documentation of this file.

// === AUDIT STATUS ===
// internal:    { status: Complete, auditors: [Khashayar], commit: }
// external_1:  { status: not started, auditors: [], commit: }
// external_2:  { status: not started, auditors: [], commit: }
// =====================
 
#pragma once
 
#include <array>
#include <span>
#include <vector>
 
#include "barretenberg/common/assert.hpp"
#include "barretenberg/common/serialize.hpp"
#include "barretenberg/polynomials/barycentric.hpp"
#include "barretenberg/polynomials/univariate_coefficient_basis.hpp"
 
namespace bb {
 
template <class Fr, size_t view_domain_end> class UnivariateView;
 
template <class Fr, size_t domain_end> class Univariate {
  public:
    static constexpr size_t LENGTH = domain_end;
    using View = UnivariateView<Fr, domain_end>;
    static constexpr size_t MONOMIAL_LENGTH = LENGTH > 1 ? 2 : 1;
    using CoefficientAccumulator = UnivariateCoefficientBasis<Fr, MONOMIAL_LENGTH, true>;
 
    using value_type = Fr; // used to get the type of the elements consistently with std::array
 
    std::array<Fr, LENGTH> evaluations;
 
    Univariate() = default;
 
    explicit Univariate(const std::array<Fr, LENGTH>& evaluations)
        : evaluations(evaluations)
    {}
    ~Univariate() = default;
    Univariate(const Univariate& other) = default;
    Univariate(Univariate&& other) noexcept = default;
    Univariate& operator=(const Univariate& other) = default;
    Univariate& operator=(Univariate&& other) noexcept = default;
 
    explicit operator UnivariateCoefficientBasis<Fr, 2, true>() const
        requires(LENGTH > 1)
    {
        static_assert(domain_end >= 2);
 
        UnivariateCoefficientBasis<Fr, 2, true> result;
        result.coefficients[0] = evaluations[0];
        result.coefficients[1] = evaluations[1] - evaluations[0];
        result.coefficients[2] = evaluations[1];
        return result;
    }
 
    // Compute Lagrange coefficients of a given linear polynomial represented in monomial basis.
    template <bool has_a0_plus_a1> Univariate(const UnivariateCoefficientBasis<Fr, 2, has_a0_plus_a1>& monomial)
    {
        // For domain_end == 1 the loop below never runs and the linear coefficient a1 would be silently
        // discarded, so this constructor cannot represent a linear polynomial on a single point.
        static_assert(domain_end >= 2);
 
        Fr to_add = monomial.coefficients[1];
        evaluations[0] = monomial.coefficients[0];
        auto prev = evaluations[0];
 
        for (size_t i = 1; i < domain_end; ++i) {
            prev = prev + to_add;
            evaluations[i] = prev;
        }
    }
 
    // Compute Lagrange coefficients of a given quadratic polynomial represented in monomial basis.
    template <bool has_a0_plus_a1> Univariate(const UnivariateCoefficientBasis<Fr, 3, has_a0_plus_a1>& monomial)
    {
        Fr to_add = monomial.coefficients[1];                                // a1 + a2
        Fr derivative = monomial.coefficients[2] + monomial.coefficients[2]; // 2a2
        evaluations[0] = monomial.coefficients[0];
        auto prev = evaluations[0];
 
        for (size_t i = 1; i < domain_end - 1; ++i) {
            prev += to_add;
            evaluations[i] = prev;
            to_add += derivative;
        }
        prev += to_add;
        evaluations[domain_end - 1] = prev;
    }
 
    // Construct constant Univariate from scalar which represents the value that all the points in the domain
    // evaluate to
    explicit Univariate(const Fr& value)
    {
        for (size_t i = 0; i < LENGTH; ++i) {
            evaluations[i] = value;
        }
    }
    // Construct Univariate from UnivariateView.
    // Lengths will match since we use `domain_end` both in the Univariate and the UnivariateView.
    explicit Univariate(const UnivariateView<Fr, domain_end>& in)
    {
        static_assert(UnivariateView<Fr, domain_end>::LENGTH == LENGTH);
        for (size_t i = 0; i < LENGTH; ++i) {
            evaluations[i] = in.evaluations[i];
        }
    }
 
    Fr& value_at(size_t i) { return evaluations[i]; }
    const Fr& value_at(size_t i) const { return evaluations[i]; }
    size_t size() { return evaluations.size(); };
 
    // Check if the univariate is identically zero
    bool is_zero() const
    {
        for (size_t i = 0; i < LENGTH; ++i) {
            // Fr::is_zero() returns bool for both scalar Fr and VectorField (the new VectorField::is_zero
            // returns true only when ALL lanes are zero — see VectorField for rationale and is_zero_mask).
            if (!evaluations[i].is_zero()) {
                return false;
            }
        }
        return true;
    }
 
    // Write the Univariate evaluations to a buffer
    [[nodiscard]] std::vector<uint8_t> to_buffer() const { return ::to_buffer(evaluations); }
 
    // Static method for creating a Univariate from a buffer
    // IMPROVEMENT: Could be made to identically match equivalent methods in e.g. field.hpp. Currently bypasses
    // unnecessary ::from_buffer call
    static Univariate serialize_from_buffer(uint8_t const* buffer)
    {
        Univariate result;
        std::read(buffer, result.evaluations);
        return result;
    }
 
    static Univariate get_random()
    {
        auto output = Univariate<Fr, domain_end>();
        for (size_t i = 0; i != LENGTH; ++i) {
            output.value_at(i) = Fr::random_element();
        }
        return output;
    };
 
    static Univariate zero() { return Univariate<Fr, domain_end>(Fr::zero()); };
 
    // Operations between Univariate and other Univariate
    bool operator==(const Univariate& other) const = default;
 
    Univariate& operator+=(const Univariate& other)
    {
        for (size_t i = 0; i < LENGTH; ++i) {
            evaluations[i] += other.evaluations[i];
        }
        return *this;
    }
    Univariate& operator-=(const Univariate& other)
    {
        for (size_t i = 0; i < LENGTH; ++i) {
            evaluations[i] -= other.evaluations[i];
        }
        return *this;
    }
    Univariate& operator*=(const Univariate& other)
    {
        for (size_t i = 0; i < LENGTH; ++i) {
            evaluations[i] *= other.evaluations[i];
        }
        return *this;
    }
    Univariate& self_sqr()
    {
        for (size_t i = 0; i < LENGTH; ++i) {
            evaluations[i].self_sqr();
        }
        return *this;
    }
    Univariate operator+(const Univariate& other) const
    {
        Univariate res(*this);
        res += other;
        return res;
    }
 
    Univariate operator-(const Univariate& other) const
    {
        Univariate res(*this);
        res -= other;
        return res;
    }
    Univariate operator-() const
    {
        Univariate res(*this);
        for (auto& eval : res.evaluations) {
            eval = -eval;
        }
        return res;
    }
 
    Univariate operator*(const Univariate& other) const
    {
        Univariate res(*this);
        res *= other;
        return res;
    }
 
    Univariate sqr() const
    {
        Univariate res(*this);
        res.self_sqr();
        return res;
    }
 
    // Operations between Univariate and scalar
    Univariate& operator+=(const Fr& scalar)
    {
        for (auto& eval : evaluations) {
            eval += scalar;
        }
        return *this;
    }
 
    Univariate& operator-=(const Fr& scalar)
    {
        for (auto& eval : evaluations) {
            eval -= scalar;
        }
        return *this;
    }
    Univariate& operator*=(const Fr& scalar)
    {
        for (auto& eval : evaluations) {
            eval *= scalar;
        }
        return *this;
    }
 
    Univariate operator+(const Fr& scalar) const
    {
        Univariate res(*this);
        res += scalar;
        return res;
    }
 
    Univariate operator-(const Fr& scalar) const
    {
        Univariate res(*this);
        res -= scalar;
        return res;
    }
 
    Univariate operator*(const Fr& scalar) const
    {
        Univariate res(*this);
        res *= scalar;
        return res;
    }
 
    // Operations between Univariate and UnivariateView
    Univariate& operator+=(const UnivariateView<Fr, domain_end>& view)
    {
        for (size_t i = 0; i < LENGTH; ++i) {
            evaluations[i] += view.evaluations[i];
        }
        return *this;
    }
 
    Univariate& operator-=(const UnivariateView<Fr, domain_end>& view)
    {
        for (size_t i = 0; i < LENGTH; ++i) {
            evaluations[i] -= view.evaluations[i];
        }
        return *this;
    }
 
    Univariate& operator*=(const UnivariateView<Fr, domain_end>& view)
    {
        for (size_t i = 0; i < LENGTH; ++i) {
            evaluations[i] *= view.evaluations[i];
        }
        return *this;
    }
 
    Univariate operator+(const UnivariateView<Fr, domain_end>& view) const
    {
        Univariate res(*this);
        res += view;
        return res;
    }
 
    Univariate operator-(const UnivariateView<Fr, domain_end>& view) const
    {
        Univariate res(*this);
        res -= view;
        return res;
    }
 
    Univariate operator*(const UnivariateView<Fr, domain_end>& view) const
    {
        Univariate res(*this);
        res *= view;
        return res;
    }
 
    // Output is immediately parsable as a list of integers by Python.
    friend std::ostream& operator<<(std::ostream& os, const Univariate& u)
    {
        os << "[";
        os << u.evaluations[0] << "," << std::endl;
        for (size_t i = 1; i < u.evaluations.size(); i++) {
            os << " " << u.evaluations[i];
            if (i + 1 < u.evaluations.size()) {
                os << "," << std::endl;
            } else {
                os << "]";
            };
        }
        return os;
    }
 
    template <size_t EXTENDED_DOMAIN_END>
    explicit operator Univariate<Fr, EXTENDED_DOMAIN_END>()
        requires(domain_end == 2)
    {
        return extend_to<EXTENDED_DOMAIN_END>();
    }
 
    template <size_t EXTENDED_DOMAIN_END> Univariate<Fr, EXTENDED_DOMAIN_END> extend_to() const
    {
        static constexpr size_t EXTENDED_LENGTH = EXTENDED_DOMAIN_END;
        using Data = BarycentricData<Fr, LENGTH, EXTENDED_LENGTH>;
        static_assert(EXTENDED_LENGTH >= LENGTH);
 
        Univariate<Fr, EXTENDED_LENGTH> result;
 
        std::copy(evaluations.begin(), evaluations.end(), result.evaluations.begin());
 
        if constexpr (LENGTH == 2) {
            // f = b x + c
            // f(0) = c
            // f(1) = b + c
            // Hence, b = f(1) - f(0)
            // f(i + 1) = f(i) + b
            Fr delta = value_at(1) - value_at(0);
            static_assert(EXTENDED_LENGTH != 0);
            for (size_t idx = domain_end - 1; idx < EXTENDED_DOMAIN_END - 1; idx++) {
                result.value_at(idx + 1) = result.value_at(idx) + delta;
            }
        } else if constexpr (LENGTH == 3) {
            static constexpr Fr inverse_two = Fr(2).invert();
            // f = a x^2 + b x + c
            // f(0) = c
            // f(1) = a + b + c
            // f(2) = 4a + 2b + c
            // f(2) + f(0) - 2f(1) = 2a
            // Hence, a = (f(2) + f(0) - 2f(1)) / 2
            // b = f(1) - a - f(0)
            // f(i+1) = f(i) + 2a * i + b + a
            // Cost note: after computing a,b, extending several points costs a few adds per point (no
            // inversions), vs the generic barycentric path.
            Fr a = (value_at(2) + value_at(0)) * inverse_two - value_at(1);
            Fr b = value_at(1) - a - value_at(0);
            Fr a2 = a + a;
            Fr a_mul = a2;
            // compute 2a * (domain_end - 1)
            for (size_t i = 0; i < domain_end - 2; i++) {
                a_mul += a2;
            }
            Fr extra = a_mul + a + b;
            for (size_t idx = domain_end - 1; idx < EXTENDED_DOMAIN_END - 1; idx++) {
                result.value_at(idx + 1) = result.value_at(idx) + extra;
                extra += a2;
            }
        } else if constexpr (LENGTH == 4) {
            static constexpr Fr inverse_six = Fr(6).invert(); // computed at compile time for efficiency
 
            // To compute a barycentric extension, we can compute the coefficients of the univariate.
            // We have the evaluation of the polynomial at the domain (which is assumed to be 0, 1, 2, 3).
            // Therefore, we have the 4 linear equations from plugging into f(x) = ax^3 + bx^2 + cx + d:
            //          a*0 + b*0 + c*0 + d = f(0)
            //          a*1 + b*1 + c*1 + d = f(1)
            //          a*2^3 + b*2^2 + c*2 + d = f(2)
            //          a*3^3 + b*3^2 + c*3 + d = f(3)
            // These equations can be rewritten as a matrix equation M * [a, b, c, d] = [f(0), f(1), f(2),
            // f(3)], where M is:
            //          0,  0,  0,  1
            //          1,  1,  1,  1
            //          2^3, 2^2, 2,  1
            //          3^3, 3^2, 3,  1
            // We can invert this matrix in order to compute a, b, c, d:
            //      -1/6,   1/2,    -1/2,   1/6
            //      1,      -5/2,   2,      -1/2
            //      -11/6,  3,      -3/2,   1/3
            //      1,      0,      0,      0
            // To compute these values, we can multiply everything by 6 and multiply by inverse_six at the
            // end for each coefficient The resulting computation here does 18 field adds, 6 subtracts, 3
            // muls to compute a, b, c, and d.
            Fr zero_times_3 = value_at(0) + value_at(0) + value_at(0);
            Fr zero_times_6 = zero_times_3 + zero_times_3;
            Fr zero_times_12 = zero_times_6 + zero_times_6;
            Fr one_times_3 = value_at(1) + value_at(1) + value_at(1);
            Fr one_times_6 = one_times_3 + one_times_3;
            Fr two_times_3 = value_at(2) + value_at(2) + value_at(2);
            Fr three_times_2 = value_at(3) + value_at(3);
            Fr three_times_3 = three_times_2 + value_at(3);
 
            Fr one_minus_two_times_3 = one_times_3 - two_times_3;
            Fr one_minus_two_times_6 = one_minus_two_times_3 + one_minus_two_times_3;
            Fr one_minus_two_times_12 = one_minus_two_times_6 + one_minus_two_times_6;
            Fr a = (one_minus_two_times_3 + value_at(3) - value_at(0)) * inverse_six; // compute a in 1 muls and 4 adds
            Fr b = (zero_times_6 - one_minus_two_times_12 - one_times_3 - three_times_3) * inverse_six;
            Fr c = (value_at(0) - zero_times_12 + one_minus_two_times_12 + one_times_6 + two_times_3 + three_times_2) *
                   inverse_six;
 
            // Then, outside of the a, b, c, d computation, we need to do some extra precomputation
            // This work is 3 field muls, 8 adds
            Fr a_plus_b = a + b;
            Fr a_plus_b_times_2 = a_plus_b + a_plus_b;
            size_t start_idx_sqr = (domain_end - 1) * (domain_end - 1);
            size_t idx_sqr_three = start_idx_sqr + start_idx_sqr + start_idx_sqr;
            Fr idx_sqr_three_times_a = Fr(idx_sqr_three) * a;
            Fr x_a_term = Fr(6 * (domain_end - 1)) * a;
            Fr three_a = a + a + a;
            Fr six_a = three_a + three_a;
 
            Fr three_a_plus_two_b = a_plus_b_times_2 + a;
            Fr linear_term = Fr(domain_end - 1) * three_a_plus_two_b + (a_plus_b + c);
            // For each new evaluation, we do only 6 field additions and 0 muls.
            for (size_t idx = domain_end - 1; idx < EXTENDED_DOMAIN_END - 1; idx++) {
                result.value_at(idx + 1) = result.value_at(idx) + idx_sqr_three_times_a + linear_term;
 
                idx_sqr_three_times_a += x_a_term + three_a;
                x_a_term += six_a;
 
                linear_term += three_a_plus_two_b;
            }
        } else {
            for (size_t k = domain_end; k != EXTENDED_DOMAIN_END; ++k) {
                result.value_at(k) = 0;
                // compute each term v_j / (d_j*(x-x_j)) of the sum
                for (size_t j = 0; j != domain_end; ++j) {
                    Fr term = value_at(j);
                    term *= Data::precomputed_denominator_inverses[LENGTH * k + j];
                    result.value_at(k) += term;
                }
                // scale the sum by the value of of B(x)
                result.value_at(k) *= Data::full_numerator_values[k];
            }
        }
        return result;
    }
 
    template <size_t INITIAL_LENGTH> void self_extend_from()
    {
        if constexpr (INITIAL_LENGTH == 2) {
            const Fr delta = value_at(1) - value_at(0);
            Fr next = value_at(1);
            for (size_t idx = 2; idx < LENGTH; idx++) {
                next += delta;
                value_at(idx) = next;
            }
        } else {
            throw_or_abort("self_extend_from called with INITIAL_LENGTH different from 2.");
        }
    }
 
    Fr evaluate(const Fr& u) const
    {
        using Data = BarycentricData<Fr, domain_end, LENGTH>;
        Fr full_numerator_value = 1;
        for (size_t i = 0; i != domain_end; ++i) {
            full_numerator_value *= u - i;
        }
 
        // 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
        std::array<Fr, LENGTH> denominator_inverses;
        for (size_t i = 0; i != LENGTH; ++i) {
            Fr inv = Data::lagrange_denominators[i];
            inv *= u - Data::big_domain[i]; // warning: need to avoid zero here
            inv = Fr(1) / inv;
            denominator_inverses[i] = inv;
        }
 
        Fr result = 0;
        // compute each term v_j / (d_j*(x-x_j)) of the sum
        for (size_t i = 0; i != domain_end; ++i) {
            Fr term = value_at(i);
            term *= denominator_inverses[i];
            result += term;
        }
        // scale the sum by the value of of B(x)
        result *= full_numerator_value;
        return result;
    };
 
    // Begin iterators
    auto begin() { return evaluations.begin(); }
    auto begin() const { return evaluations.begin(); }
    // End iterators
    auto end() { return evaluations.end(); }
    auto end() const { return evaluations.end(); }
};
 
template <typename B, class Fr, size_t domain_end> inline void read(B& it, Univariate<Fr, domain_end>& univariate)
{
    using serialize::read;
    read(it, univariate.evaluations);
}
 
template <typename B, class Fr, size_t domain_end>
inline void write(B& it, Univariate<Fr, domain_end> const& univariate)
{
    using serialize::write;
    write(it, univariate.evaluations);
}
 
template <class Fr, size_t domain_end>
Univariate<Fr, domain_end> operator+(const Fr& ff, const Univariate<Fr, domain_end>& uv)
{
    return uv + ff;
}
 
template <class Fr, size_t domain_end>
Univariate<Fr, domain_end> operator-(const Fr& ff, const Univariate<Fr, domain_end>& uv)
{
    return -uv + ff;
}
 
template <class Fr, size_t domain_end>
Univariate<Fr, domain_end> operator*(const Fr& ff, const Univariate<Fr, domain_end>& uv)
{
    return uv * ff;
}
 
template <class Fr, size_t domain_end> class UnivariateView {
  public:
    static constexpr size_t LENGTH = domain_end;
    std::span<const Fr, LENGTH> evaluations;
    static constexpr size_t MONOMIAL_LENGTH = LENGTH > 1 ? 2 : 1;
    using CoefficientAccumulator = UnivariateCoefficientBasis<Fr, MONOMIAL_LENGTH, true>;
 
    UnivariateView() = default;
 
    bool operator==(const UnivariateView& other) const
    {
        for (size_t i = 0; i < LENGTH; ++i) {
            if (evaluations[i] != other.evaluations[i]) {
                return false;
            }
        }
        return true;
    };
 
    const Fr& value_at(size_t i) const { return evaluations[i]; };
 
    template <size_t full_domain_end>
    explicit UnivariateView(const Univariate<Fr, full_domain_end>& univariate_in)
        : evaluations(std::span<const Fr>(univariate_in.evaluations.data(), LENGTH))
    {
        // Viewing more evaluations than the source holds would span out-of-bounds memory.
        static_assert(LENGTH <= full_domain_end);
    };
 
    explicit operator UnivariateCoefficientBasis<Fr, 2, true>() const
        requires(LENGTH > 1)
    {
        static_assert(domain_end >= 2);
 
        UnivariateCoefficientBasis<Fr, 2, true> result;
 
        result.coefficients[0] = evaluations[0];
        result.coefficients[1] = evaluations[1] - evaluations[0];
        result.coefficients[2] = evaluations[1];
        return result;
    }
 
    Univariate<Fr, domain_end> operator+(const UnivariateView& other) const
    {
        Univariate<Fr, domain_end> res(*this);
        res += other;
        return res;
    }
 
    Univariate<Fr, domain_end> operator-(const UnivariateView& other) const
    {
        Univariate<Fr, domain_end> res(*this);
        res -= other;
        return res;
    }
 
    Univariate<Fr, domain_end> operator-() const
    {
        Univariate<Fr, domain_end> res(*this);
        for (auto& eval : res.evaluations) {
            eval = -eval;
        }
        return res;
    }
 
    Univariate<Fr, domain_end> operator*(const UnivariateView& other) const
    {
        Univariate<Fr, domain_end> res(*this);
        res *= other;
        return res;
    }
    Univariate<Fr, domain_end> sqr() const
    {
        Univariate<Fr, domain_end> res(*this);
        res = res.sqr();
        return res;
    }
 
    Univariate<Fr, domain_end> operator*(const Univariate<Fr, domain_end>& other) const
    {
        Univariate<Fr, domain_end> res(*this);
        res *= other;
        return res;
    }
 
    Univariate<Fr, domain_end> operator+(const Univariate<Fr, domain_end>& other) const
    {
        Univariate<Fr, domain_end> res(*this);
        res += other;
        return res;
    }
 
    Univariate<Fr, domain_end> operator+(const Fr& other) const
    {
        Univariate<Fr, domain_end> res(*this);
        res += other;
        return res;
    }
 
    Univariate<Fr, domain_end> operator-(const Fr& other) const
    {
        Univariate<Fr, domain_end> res(*this);
        res -= other;
        return res;
    }
 
    Univariate<Fr, domain_end> operator*(const Fr& other) const
    {
        Univariate<Fr, domain_end> res(*this);
        res *= other;
        return res;
    }
 
    Univariate<Fr, domain_end> operator-(const Univariate<Fr, domain_end>& other) const
    {
        Univariate<Fr, domain_end> res(*this);
        res -= other;
        return res;
    }
 
    // Output is immediately parsable as a list of integers by Python.
    friend std::ostream& operator<<(std::ostream& os, const UnivariateView& u)
    {
        os << "[";
        os << u.evaluations[0] << "," << std::endl;
        for (size_t i = 1; i < u.evaluations.size(); i++) {
            os << " " << u.evaluations[i];
            if (i + 1 < u.evaluations.size()) {
                os << "," << std::endl;
            } else {
                os << "]";
            };
        }
        return os;
    }
};
 
template <class Fr, size_t domain_end>
Univariate<Fr, domain_end> operator+(const Fr& ff, const UnivariateView<Fr, domain_end>& uv)
{
    return uv + ff;
}
 
template <class Fr, size_t domain_end>
Univariate<Fr, domain_end> operator-(const Fr& ff, const UnivariateView<Fr, domain_end>& uv)
{
    return -uv + ff;
}
 
template <class Fr, size_t domain_end>
Univariate<Fr, domain_end> operator*(const Fr& ff, const UnivariateView<Fr, domain_end>& uv)
{
    return uv * ff;
}
 
template <typename T, typename U, std::size_t N, std::size_t... Is>
std::array<T, sizeof...(Is)> array_to_array_aux(const std::array<U, N>& elements, std::index_sequence<Is...>)
{
    return { { T{ elements[Is] }... } };
};
 
template <typename T, typename U, std::size_t N> std::array<T, N> array_to_array(const std::array<U, N>& elements)
{
    // Calls the aux method that uses the index sequence to unpack all values in `elements`
    return array_to_array_aux<T, U, N>(elements, std::make_index_sequence<N>());
};
 
} // namespace bb
 
namespace std {
 
template <typename T, size_t N> struct tuple_size<bb::Univariate<T, N>> : std::integral_constant<std::size_t, N> {};
 
} // namespace std