biggroup_impl.hpp

Go to the documentation of this file.

// === AUDIT STATUS ===
// internal:    { status: not started, auditors: [], date: YYYY-MM-DD }
// external_1:  { status: not started, auditors: [], date: YYYY-MM-DD }
// external_2:  { status: not started, auditors: [], date: YYYY-MM-DD }
// =====================
 
#pragma once
 
#include "../circuit_builders/circuit_builders.hpp"
#include "../plookup/plookup.hpp"
#include "barretenberg/common/assert.hpp"
#include "barretenberg/ecc/groups/precomputed_generators.hpp"
#include "barretenberg/numeric/general/general.hpp"
#include "barretenberg/stdlib/primitives/biggroup/biggroup.hpp"
#include "barretenberg/transcript/origin_tag.hpp"
 
namespace bb::stdlib::element_default {
 
template <typename C, class Fq, class Fr, class G>
element<C, Fq, Fr, G>::element()
    : _x()
    , _y()
    , _is_infinity()
{}
 
template <typename C, class Fq, class Fr, class G>
element<C, Fq, Fr, G>::element(const typename G::affine_element& input)
    : _x(nullptr, input.x)
    , _y(nullptr, input.y)
    , _is_infinity(nullptr, input.is_point_at_infinity())
{}
 
template <typename C, class Fq, class Fr, class G>
element<C, Fq, Fr, G>::element(const Fq& x_in, const Fq& y_in, const bool assert_on_curve)
    : _x(x_in)
    , _y(y_in)
    , _is_infinity(_x.get_context() ? _x.get_context() : _y.get_context(), false)
{
    if (assert_on_curve) {
        validate_on_curve();
    }
}
 
template <typename C, class Fq, class Fr, class G>
element<C, Fq, Fr, G>::element(const Fq& x_in, const Fq& y_in, const bool_ct& is_infinity, const bool assert_on_curve)
    : _x(x_in)
    , _y(y_in)
    , _is_infinity(is_infinity)
{
    if (assert_on_curve) {
        validate_on_curve();
    }
}
 
template <typename C, class Fq, class Fr, class G>
element<C, Fq, Fr, G>::element(const element& other)
    : _x(other._x)
    , _y(other._y)
    , _is_infinity(other.is_point_at_infinity())
{}
 
template <typename C, class Fq, class Fr, class G>
element<C, Fq, Fr, G>::element(element&& other) noexcept
    : _x(other._x)
    , _y(other._y)
    , _is_infinity(other.is_point_at_infinity())
{}
 
template <typename C, class Fq, class Fr, class G>
element<C, Fq, Fr, G>& element<C, Fq, Fr, G>::operator=(const element& other)
{
    if (&other == this) {
        return *this;
    }
    _x = other._x;
    _y = other._y;
    _is_infinity = other.is_point_at_infinity();
    return *this;
}
 
template <typename C, class Fq, class Fr, class G>
element<C, Fq, Fr, G>& element<C, Fq, Fr, G>::operator=(element&& other) noexcept
{
    if (&other == this) {
        return *this;
    }
    _x = other._x;
    _y = other._y;
    _is_infinity = other.is_point_at_infinity();
    return *this;
}
 
template <typename C, class Fq, class Fr, class G>
element<C, Fq, Fr, G> element<C, Fq, Fr, G>::operator+(const element& other) const
{
    // Adding in `x_coordinates_match` ensures that lambda will always be well-formed
    // Our curve has the form y² = x³ + ax + b (or y² = x³ + b when a = 0).
    // If (x₁, y₁), (x₂, y₂) have x₁ == x₂, the generic formula for lambda has a division by 0.
    // Then y₁ == y₂ (i.e. we are doubling) or y₂ == -y₁ (the sum is infinity).
    // These cases have special addition formulae. The following booleans allow us to handle these cases uniformly.
    const bool_ct x_coordinates_match = other._x == _x;
    const bool_ct y_coordinates_match = (_y == other._y);
    const bool_ct infinity_predicate = (x_coordinates_match && !y_coordinates_match);
    const bool_ct double_predicate = (x_coordinates_match && y_coordinates_match);
    const bool_ct lhs_infinity = is_point_at_infinity();
    const bool_ct rhs_infinity = other.is_point_at_infinity();
    const bool_ct has_infinity_input = lhs_infinity || rhs_infinity;
 
    // NOTE: For valid points on the curve, specifically for bn254 or secp256k1 or secp256r1, y = 0 cannot occur.
    // For points not on the curve, having y = 0 will lead to a failure while performing the division below.
    // We could enforce in circuit that y = 0 results in point at infinity, or that y != 0 always.
    // However, this would be an unnecessary constraint for valid points on the curve.
    // So we perform a native check here to catch any accidental misuse of this function.
    const typename G::Fq y_value = uint256_t(_y.get_value());
    BB_ASSERT_EQ((y_value == 0), false, "Attempting to add a point with y = 0, not allowed.");
 
    // Compute the gradient λ. If we add, λ = (y₂ - y₁)/(x₂ - x₁)
    // For doubling: λ = (3x₁² + a)/(2y₁) if curve has 'a', else λ = 3x₁²/(2y₁)
    const Fq add_lambda_numerator = other._y - _y;
    const Fq xx = _x * _x;
    Fq dbl_lambda_numerator = xx + xx + xx; // 3x²
    if constexpr (G::has_a) {
        // Curve equation: y² = x³ + ax + b
        // Doubling formula numerator: 3x² + a
        const Fq a(get_context(), uint256_t(G::curve_a));
        dbl_lambda_numerator = dbl_lambda_numerator + a;
    }
    const Fq lambda_numerator = Fq::conditional_assign(double_predicate, dbl_lambda_numerator, add_lambda_numerator);
 
    const Fq add_lambda_denominator = other._x - _x;
    const Fq dbl_lambda_denominator = _y + _y;
    Fq lambda_denominator = Fq::conditional_assign(double_predicate, dbl_lambda_denominator, add_lambda_denominator);
 
    // If either input is a point at infinity, or if the result would be infinity, set lambda_denominator to 1
    // to prevent division by zero. Cases where result is infinity: x₁ == x₂ but y₁ != y₂ (points are inverses)
    const bool_ct safe_denominator_needed = has_infinity_input || infinity_predicate;
    lambda_denominator = Fq::conditional_assign(safe_denominator_needed, Fq(1), lambda_denominator);
    const Fq lambda = Fq::div_without_denominator_check({ lambda_numerator }, lambda_denominator);
 
    // Compute resulting point coordinates: x₃ = λ² - x₁ - x₂, y₃ = λ(x₁ - x₃) - y₁
    const Fq x3 = lambda.sqradd({ -other._x, -_x });
    const Fq y3 = lambda.madd(_x - x3, { -_y });
    element result(x3, y3, /*assert_on_curve=*/false);
 
    // if lhs infinity, return rhs
    result._x = Fq::conditional_assign(lhs_infinity, other._x, result._x);
    result._y = Fq::conditional_assign(lhs_infinity, other._y, result._y);
    // if rhs infinity, return lhs
    result._x = Fq::conditional_assign(rhs_infinity, _x, result._x);
    result._y = Fq::conditional_assign(rhs_infinity, _y, result._y);
 
    // Determine if result is point at infinity:
    // - If x₁ == x₂ and y₁ == -y₂ (i.e., points are inverses), result is ∞
    // - If both inputs are ∞, result is ∞
    bool_ct result_is_infinity = (infinity_predicate && !has_infinity_input) || (lhs_infinity && rhs_infinity);
    result.set_point_at_infinity(result_is_infinity, /* add_to_used_witnesses */ true);
 
    result.set_origin_tag(OriginTag(get_origin_tag(), other.get_origin_tag()));
    return result;
}
 
template <typename C, class Fq, class Fr, class G>
element<C, Fq, Fr, G> element<C, Fq, Fr, G>::get_standard_form() const
{
 
    const bool_ct is_infinity = is_point_at_infinity();
    element result(*this);
    const Fq zero = Fq::zero();
    result._x = Fq::conditional_assign(is_infinity, zero, this->_x);
    result._y = Fq::conditional_assign(is_infinity, zero, this->_y);
    return result;
}
 
template <typename C, class Fq, class Fr, class G>
element<C, Fq, Fr, G> element<C, Fq, Fr, G>::operator-(const element& other) const
{
    // Adding in `x_coordinates_match` ensures that lambda will always be well-formed
    // Our curve has the form y² = x³ + ax + b (or y² = x³ + b when a = 0).
    // If (x₁, y₁), (x₂, y₂) have x₁ == x₂, the generic formula for lambda has a division by 0.
    // For subtraction P₁ - P₂ = P₁ + (-P₂), where -P₂ = (x₂, -y₂):
    // - If y₁ == -y₂ (i.e., x₁ == x₂ and y₁ == -y₂), this becomes doubling: P₁ + P₁
    // - If y₁ == y₂ (i.e., x₁ == x₂ and y₁ == y₂), result is infinity: P₁ - P₁ = ∞
    // These cases have special addition formulae. The following booleans allow us to handle these cases uniformly.
    const bool_ct x_coordinates_match = other._x == _x;
    const bool_ct y_coordinates_match = (_y == other._y);
    const bool_ct infinity_predicate = (x_coordinates_match && y_coordinates_match);
    const bool_ct double_predicate = (x_coordinates_match && !y_coordinates_match);
    const bool_ct lhs_infinity = is_point_at_infinity();
    const bool_ct rhs_infinity = other.is_point_at_infinity();
    const bool_ct has_infinity_input = lhs_infinity || rhs_infinity;
 
    // Compute the gradient λ. For subtraction, λ = (-y₂ - y₁)/(x₂ - x₁)
    // For doubling: λ = (3x₁² + a)/(2y₁) if curve has 'a', else λ = 3x₁²/(2y₁)
    const Fq add_lambda_numerator = -other._y - _y;
    const Fq xx = _x * _x;
    Fq dbl_lambda_numerator = xx + xx + xx; // 3x²
    if constexpr (G::has_a) {
        // Curve equation: y² = x³ + ax + b
        // Doubling formula numerator: 3x² + a
        const Fq a(get_context(), uint256_t(G::curve_a));
        dbl_lambda_numerator = dbl_lambda_numerator + a;
    }
    const Fq lambda_numerator = Fq::conditional_assign(double_predicate, dbl_lambda_numerator, add_lambda_numerator);
 
    const Fq add_lambda_denominator = other._x - _x;
    const Fq dbl_lambda_denominator = _y + _y;
    Fq lambda_denominator = Fq::conditional_assign(double_predicate, dbl_lambda_denominator, add_lambda_denominator);
 
    // If either input is a point at infinity, or if the result would be infinity (x₁ == x₂ and y₁ == y₂),
    // set lambda_denominator to 1 to prevent division by zero. The lambda value won't be used in these cases.
    lambda_denominator = Fq::conditional_assign(has_infinity_input || infinity_predicate, Fq(1), lambda_denominator);
    const Fq lambda = Fq::div_without_denominator_check({ lambda_numerator }, lambda_denominator);
 
    // Compute resulting point coordinates: x₃ = λ² - x₁ - x₂, y₃ = λ(x₁ - x₃) - y₁
    const Fq x3 = lambda.sqradd({ -other._x, -_x });
    const Fq y3 = lambda.madd(_x - x3, { -_y });
    element result(x3, y3, /*assert_on_curve=*/false);
 
    // if lhs infinity, return -rhs (negated rhs point)
    result._x = Fq::conditional_assign(lhs_infinity, other._x, result._x);
    result._y = Fq::conditional_assign(lhs_infinity, -other._y, result._y);
    // if rhs infinity, return lhs
    result._x = Fq::conditional_assign(rhs_infinity, _x, result._x);
    result._y = Fq::conditional_assign(rhs_infinity, _y, result._y);
 
    // Determine if result is point at infinity:
    // - If x₁ == x₂ and y₁ == y₂ (i.e., P₁ - P₁), result is ∞
    // - If both inputs are ∞, result is ∞
    bool_ct result_is_infinity = (infinity_predicate && !has_infinity_input) || (lhs_infinity && rhs_infinity);
 
    result.set_point_at_infinity(result_is_infinity, /* add_to_used_witnesses */ true);
    result.set_origin_tag(OriginTag(get_origin_tag(), other.get_origin_tag()));
    return result;
}
 
template <typename C, class Fq, class Fr, class G>
element<C, Fq, Fr, G> element<C, Fq, Fr, G>::checked_unconditional_add(const element& other) const
{
    other._x.assert_is_not_equal(_x);
    const Fq lambda = Fq::div_without_denominator_check({ other._y, -_y }, (other._x - _x));
    const Fq x3 = lambda.sqradd({ -other._x, -_x });
    const Fq y3 = lambda.madd(_x - x3, { -_y });
    return element(x3, y3, /*assert_on_curve=*/false);
}
 
template <typename C, class Fq, class Fr, class G>
element<C, Fq, Fr, G> element<C, Fq, Fr, G>::checked_unconditional_subtract(const element& other) const
{
 
    other._x.assert_is_not_equal(_x);
    const Fq lambda = Fq::div_without_denominator_check({ other._y, _y }, (other._x - _x));
    const Fq x_3 = lambda.sqradd({ -other._x, -_x });
    const Fq y_3 = lambda.madd(x_3 - _x, { -_y });
 
    return element(x_3, y_3, /*assert_on_curve=*/false);
}
 
template <typename C, class Fq, class Fr, class G>
std::array<element<C, Fq, Fr, G>, 2> element<C, Fq, Fr, G>::checked_unconditional_add_sub(const element& other) const
{
    // validate we can use incomplete addition formulae
    other._x.assert_is_not_equal(_x);
 
    const Fq denominator = other._x - _x;
    const Fq x2x1 = -(other._x + _x);
 
    const Fq lambda1 = Fq::div_without_denominator_check({ other._y, -_y }, denominator);
    const Fq x_3 = lambda1.sqradd({ x2x1 });
    const Fq y_3 = lambda1.madd(_x - x_3, { -_y });
    const Fq lambda2 = Fq::div_without_denominator_check({ -other._y, -_y }, denominator);
    const Fq x_4 = lambda2.sqradd({ x2x1 });
    const Fq y_4 = lambda2.madd(_x - x_4, { -_y });
 
    return { element(x_3, y_3, /*assert_on_curve=*/false), element(x_4, y_4, /*assert_on_curve=*/false) };
}
 
template <typename C, class Fq, class Fr, class G> element<C, Fq, Fr, G> element<C, Fq, Fr, G>::dbl() const
{
    // NOTE: For valid points on the curve, specifically for bn254 or secp256k1 or secp256r1, y = 0 cannot occur.
    // For points not on the curve, having y = 0 will lead to a failure while performing the division below.
    // We could enforce in circuit that y = 0 results in point at infinity, or that y != 0 always.
    // However, this would be an unnecessary constraint for valid points on the curve.
    // So we perform a native check here to catch any accidental misuse of this function.
    const typename G::Fq y_value = uint256_t(_y.get_value());
    BB_ASSERT_EQ((y_value == 0), false, "Attempting to dbl a point with y = 0, not allowed.");
 
    Fq two_x = _x + _x;
    if constexpr (G::has_a) {
        // Curve equation: y² = x³ + ax + b
        Fq a(get_context(), uint256_t(G::curve_a));
 
        // Compute neg_lambda = -λ = -(3x² + a) / (2y)
        // msub_div computes: -(Σᵢ aᵢ·bᵢ + Σⱼ cⱼ) / d = -(x·(3x) + a) / (2y) = -(3x² + a) / (2y)
        Fq neg_lambda = Fq::msub_div({ _x }, { (two_x + _x) }, (_y + _y), { a }, /*enable_divisor_nz_check*/ false);
 
        // Compute x₃ = λ² - 2x
        // Since neg_lambda = -λ, we have: (-λ)² - 2x = λ² - 2x
        Fq x_3 = neg_lambda.sqradd({ -(two_x) });
 
        // Compute y₃ = λ(x - x₃) - y
        // Using neg_lambda = -λ: (-λ)(x₃ - x) + (-y) = -λ(x₃ - x) - y = λ(x - x₃) - y
        Fq y_3 = neg_lambda.madd(x_3 - _x, { -_y });
 
        element result(x_3, y_3, /*assert_on_curve=*/false);
        result.set_point_at_infinity(is_point_at_infinity(), /* add_to_used_witnesses */ true);
        return result;
    }
 
    // Curve equation when a = 0: y² = x³ + b
    // Compute neg_lambda = -λ = -3x² / (2y)
    // msub_div computes: -(Σᵢ aᵢ·bᵢ) / d = -(x·(3x)) / (2y) = -3x² / (2y)
    Fq neg_lambda = Fq::msub_div({ _x }, { (two_x + _x) }, (_y + _y), {}, /*enable_divisor_nz_check*/ false);
 
    // Compute x₃ = λ² - 2x
    // Since neg_lambda = -λ, we have: (-λ)² - 2x = λ² - 2x
    Fq x_3 = neg_lambda.sqradd({ -(two_x) });
 
    // Compute y₃ = λ(x - x₃) - y
    // Using neg_lambda = -λ: (-λ)(x₃ - x) + (-y) = -λ(x₃ - x) - y = λ(x - x₃) - y
    Fq y_3 = neg_lambda.madd(x_3 - _x, { -_y });
 
    element result = element(x_3, y_3, /*assert_on_curve=*/false);
    result.set_point_at_infinity(is_point_at_infinity(), /* add_to_used_witnesses */ true);
    return result;
}
 
template <typename C, class Fq, class Fr, class G>
typename element<C, Fq, Fr, G>::chain_add_accumulator element<C, Fq, Fr, G>::chain_add_start(const element& p1,
                                                                                             const element& p2)
{
    chain_add_accumulator output;
    output.x1_prev = p1._x;
    output.y1_prev = p1._y;
 
    // Require x₁ ≠ x₂ for incomplete addition formula
    p1._x.assert_is_not_equal(p2._x);
 
    // Compute λ = (y₂ - y₁)/(x₂ - x₁)
    const Fq lambda = Fq::div_without_denominator_check({ p2._y, -p1._y }, (p2._x - p1._x));
 
    // Compute x₃ = λ² - x₁ - x₂
    const Fq x3 = lambda.sqradd({ -p2._x, -p1._x });
    output.x3_prev = x3;
    output.lambda_prev = lambda;
    return output;
}
 
template <typename C, class Fq, class Fr, class G>
typename element<C, Fq, Fr, G>::chain_add_accumulator element<C, Fq, Fr, G>::chain_add(const element& p1,
                                                                                       const chain_add_accumulator& acc)
{
    // If accumulator has a full y-coordinate, use chain_add_start instead
    if (acc.is_full_element) {
        return chain_add_start(p1, element(acc.x3_prev, acc.y3_prev, /*assert_on_curve=*/false));
    }
 
    // Require x₁ ≠ x₂ for incomplete addition formula
    p1._x.assert_is_not_equal(acc.x3_prev);
 
    // Compute λ = (y₂ - y₁)/(x₂ - x₁), but we don't have y₂!
    // We know that y₂ = lambda_prev * (x1_prev - x₂) - y1_prev from the previous addition.
    //
    // Derivation:
    //   λ(x₂ - x₁) = y₂ - y₁
    //   λ(x₂ - x₁) = lambda_prev * (x1_prev - x₂) - y1_prev - y₁
    //   λ(x₂ - x₁) = -lambda_prev * (x₂ - x1_prev) - y1_prev - y₁
    //   λ = -(lambda_prev * (x₂ - x1_prev) + y1_prev + y₁) / (x₂ - x₁)
    //
    auto& x2 = acc.x3_prev;
    const auto lambda = Fq::msub_div({ acc.lambda_prev },
                                     { (x2 - acc.x1_prev) },
                                     (x2 - p1._x),
                                     { acc.y1_prev, p1._y },
                                     /*enable_divisor_nz_check*/ false); // Divisor is non-zero as x₂ ≠ x₁ is enforced
 
    // Compute x₃ = λ² - x₂ - x₁
    const auto x3 = lambda.sqradd({ -x2, -p1._x });
 
    // Update the accumulator
    chain_add_accumulator output;
    output.x3_prev = x3;
    output.x1_prev = p1._x;
    output.y1_prev = p1._y;
    output.lambda_prev = lambda;
 
    return output;
}
 
template <typename C, class Fq, class Fr, class G>
element<C, Fq, Fr, G> element<C, Fq, Fr, G>::chain_add_end(const chain_add_accumulator& acc)
{
    // If accumulator already has a full y-coordinate, return it directly
    if (acc.is_full_element) {
        return element(acc.x3_prev, acc.y3_prev, /*assert_on_curve=*/false);
    }
 
    // Compute y₃ = λ(x₁ - x₃) - y₁
    // where λ, x₁, y₁ are from the previous addition stored in the accumulator
    auto& x3 = acc.x3_prev;
    auto& lambda = acc.lambda_prev;
 
    Fq y3 = lambda.madd((acc.x1_prev - x3), { -acc.y1_prev });
    return element(x3, y3, /*assert_on_curve=*/false);
}
 
template <typename C, class Fq, class Fr, class G>
element<C, Fq, Fr, G> element<C, Fq, Fr, G>::multiple_montgomery_ladder(
    const std::vector<chain_add_accumulator>& add) const
{
    struct composite_y {
        std::vector<Fq> mul_left;
        std::vector<Fq> mul_right;
        std::vector<Fq> add;
        bool is_negative = false;
    };
 
    // Handle edge case of empty input
    if (add.empty()) {
        return *this;
    }
 
    // Let A = (x, y) and P = (x₁, y₁)
    // For the first point P, we want to compute: (2A + P) = (A + P) + A
    // We first need to check if x ≠ x₁.
    x().assert_is_not_equal(add[0].x3_prev);
 
    // Compute λ₁ for computing the first addition: (A + P)
    Fq lambda1;
    if (!add[0].is_full_element) {
        // Case 1: P is an accumulator (i.e., it lacks a y-coordinate)
        //         λ₁ = (y - y₁) / (x - x₁)
        //            = -(y₁ - y) / (x - x₁)
        //            = -(λ₁_ₚᵣₑᵥ * (x₁_ₚᵣₑᵥ - x₁) - y₁_ₚᵣₑᵥ - y) / (x - x₁)
        //
        // NOTE: msub_div computes -(∑ᵢ aᵢ * bᵢ + ∑ⱼcⱼ) / d
        lambda1 = Fq::msub_div({ add[0].lambda_prev },              // numerator left multiplicands: λ₁_ₚᵣₑᵥ
                               { add[0].x1_prev - add[0].x3_prev }, // numerator right multiplicands: (x₁_ₚᵣₑᵥ - x₁)
                               (x() - add[0].x3_prev),              // denominator: (x - x₁)
                               { -add[0].y1_prev, -y() },           // numerator additions: -y₁_ₚᵣₑᵥ - y
                               /*enable_divisor_nz_check*/ false);  // divisor check is not needed as x ≠ x₁ is enforced
    } else {
        // Case 2: P is a full element (i.e., it has a y-coordinate)
        //         λ₁ = (y - y₁) / (x - x₁)
        //
        lambda1 = Fq::div_without_denominator_check({ y() - add[0].y3_prev }, (x() - add[0].x3_prev));
    }
 
    // Using λ₁, compute x₃ for (A + P):
    // x₃ = λ₁.λ₁ - x₁ - x
    Fq x_3 = lambda1.madd(lambda1, { -add[0].x3_prev, -x() });
 
    // Compute λ₂ for the addition (A + P) + A:
    // λ₂ = (y - y₃) / (x - x₃)
    //    = (y - (λ₁ * (x - x₃) - y)) / (x - x₃)    (substituting y₃)
    //    = (2y) / (x - x₃) - λ₁
    //
    x().assert_is_not_equal(x_3);
    Fq lambda2 = Fq::div_without_denominator_check({ y() + y() }, (x() - x_3)) - lambda1;
 
    // Using λ₂, compute x₄ for the final result:
    // x₄ = λ₂.λ₂ - x₃ - x
    Fq x_4 = lambda2.sqradd({ -x_3, -x() });
 
    // Compute y₄ for the final result:
    // y₄ = λ₂ * (x - x₄) - y
    //
    // However, we don't actually compute y₄ here. Instead, we build a "composite" y value that contains
    // the components needed to compute y₄ later. This is done to avoid the explicit multiplication here.
    //
    // We store the result as either y₄ or -y₄, depending on whether the number of points added
    // is even or odd. This sign adjustment simplifies the handling of subsequent additions in the loop below.
    // +y₄ = λ₂ * (x - x₄) - y
    // -y₄ = λ₂ * (x₄ - x) + y
    const bool num_points_even = ((add.size() & 1ULL) == 0);
    composite_y previous_y;
    previous_y.add.emplace_back(num_points_even ? y() : -y());
    previous_y.mul_left.emplace_back(lambda2);
    previous_y.mul_right.emplace_back(num_points_even ? x_4 - x() : x() - x_4);
    previous_y.is_negative = num_points_even;
 
    // Handle remaining iterations (i > 0) in a loop
    Fq previous_x = x_4;
    for (size_t i = 1; i < add.size(); ++i) {
        // Let x = previous_x, y = previous_y
        // Let P = (xᵢ, yᵢ) be the next point to add (represented by add[i])
        // Ensure x-coordinates are distinct: x ≠ xᵢ
        previous_x.assert_is_not_equal(add[i].x3_prev);
 
        // Determine sign adjustment based on previous y's sign
        // If the previous y was positive, we need to negate the y-component from add[i]
        const bool negate_add_y = !previous_y.is_negative;
 
        // Build λ₁ numerator components from previous composite y and current accumulator
        std::vector<Fq> lambda1_left = previous_y.mul_left;
        std::vector<Fq> lambda1_right = previous_y.mul_right;
        std::vector<Fq> lambda1_add = previous_y.add;
 
        if (!add[i].is_full_element) {
            // Case 1: add[i] is an accumulator (lacks y-coordinate)
            //         λ₁ = (y - yᵢ) / (x - xᵢ)
            //            = -(yᵢ - y) / (x - xᵢ)
            //            = -(λᵢ_ₚᵣₑᵥ * (xᵢ_ₚᵣₑᵥ - xᵢ) - yᵢ_ₚᵣₑᵥ - y) / (x - xᵢ)
            //
            // If (previous) y is stored as positive, we compute λ₁ as:
            //         λ₁ = -(λᵢ_ₚᵣₑᵥ * (xᵢ - xᵢ_ₚᵣₑᵥ) + yᵢ_ₚᵣₑᵥ + y) / (xᵢ - x)
            //
            lambda1_left.emplace_back(add[i].lambda_prev);
            lambda1_right.emplace_back(negate_add_y ? add[i].x3_prev - add[i].x1_prev
                                                    : add[i].x1_prev - add[i].x3_prev);
            lambda1_add.emplace_back(negate_add_y ? add[i].y1_prev : -add[i].y1_prev);
        } else {
            // Case 2: add[i] is a full element (has y-coordinate)
            //         λ₁ = (yᵢ - y) / (xᵢ - x)
            //
            // If previous y is positive, we compute λ₁ as:
            //         λ₁ = -(y - yᵢ) / (xᵢ - x)
            //
            lambda1_add.emplace_back(negate_add_y ? -add[i].y3_prev : add[i].y3_prev);
        }
 
        // Compute λ₁
        Fq denominator = negate_add_y ? add[i].x3_prev - previous_x : previous_x - add[i].x3_prev;
        Fq lambda1 =
            Fq::msub_div(lambda1_left, lambda1_right, denominator, lambda1_add, /*enable_divisor_nz_check*/ false);
 
        // Using λ₁, compute x₃ for (previous + P):
        // x₃ = λ₁.λ₁ - xᵢ - x
        // y₃ = λ₁ * (x - x₃) - y (we don't compute this explicitly)
        Fq x_3 = lambda1.madd(lambda1, { -add[i].x3_prev, -previous_x });
 
        // Compute λ₂ using previous composite y
        // λ₂ = (y - y₃) / (x - x₃)
        //    = (y - (λ₁ * (x - x₃) - y)) / (x - x₃)    (substituting y₃)
        //    = (2y) / (x - x₃) - λ₁
        //    = -2(y / (x₃ - x)) - λ₁
        //
        previous_x.assert_is_not_equal(x_3);
        Fq l2_denominator = previous_y.is_negative ? previous_x - x_3 : x_3 - previous_x;
        Fq partial_lambda2 = Fq::msub_div(previous_y.mul_left,
                                          previous_y.mul_right,
                                          l2_denominator,
                                          previous_y.add,
                                          /*enable_divisor_nz_check*/ false);
        partial_lambda2 = partial_lambda2 + partial_lambda2;
        lambda2 = partial_lambda2 - lambda1;
 
        // Using λ₂, compute x₄ for the final result of this iteration:
        // x₄ = λ₂.λ₂ - x₃ - x
        x_4 = lambda2.sqradd({ -x_3, -previous_x });
 
        // Build composite y for this iteration
        // y₄ = λ₂ * (x - x₄) - y
        // However, we don't actually compute y₄ explicitly, we rather store components to compute it later.
        // We store the result as either y₄ or -y₄, depending on the sign of previous_y.
        // +y₄ = λ₂ * (x - x₄) - y
        // -y₄ = λ₂ * (x₄ - x) + y
        composite_y y_4;
        y_4.is_negative = !previous_y.is_negative;
        y_4.mul_left.emplace_back(lambda2);
        y_4.mul_right.emplace_back(previous_y.is_negative ? previous_x - x_4 : x_4 - previous_x);
 
        // Append terms from previous_y to y_4. We want to make sure the terms above are added into the start
        // of y_4. This is to ensure they are cached correctly when
        // `builder::evaluate_partial_non_native_field_multiplication` is called. (the 1st mul_left, mul_right elements
        // will trigger builder::evaluate_non_native_field_multiplication
        //  when Fq::mult_madd is called - this term cannot be cached so we want to make sure it is unique)
        std::copy(previous_y.mul_left.begin(), previous_y.mul_left.end(), std::back_inserter(y_4.mul_left));
        std::copy(previous_y.mul_right.begin(), previous_y.mul_right.end(), std::back_inserter(y_4.mul_right));
        std::copy(previous_y.add.begin(), previous_y.add.end(), std::back_inserter(y_4.add));
 
        previous_x = x_4;
        previous_y = y_4;
    }
 
    Fq x_out = previous_x;
 
    BB_ASSERT(!previous_y.is_negative);
 
    Fq y_out = Fq::mult_madd(previous_y.mul_left, previous_y.mul_right, previous_y.add);
    return element(x_out, y_out, /*assert_on_curve=*/false);
}
 
template <typename C, class Fq, class Fr, class G>
std::pair<element<C, Fq, Fr, G>, element<C, Fq, Fr, G>> element<C, Fq, Fr, G>::compute_offset_generators(
    const size_t num_rounds)
{
    constexpr typename G::affine_element offset_generator =
        get_precomputed_generators<G, "biggroup offset generator", 1>()[0];
 
    const uint256_t offset_multiplier = uint256_t(1) << uint256_t(num_rounds - 1);
 
    const typename G::affine_element offset_generator_end = typename G::element(offset_generator) * offset_multiplier;
    return std::make_pair<element, element>(offset_generator, offset_generator_end);
}
 
template <typename C, class Fq, class Fr, class G>
element<C, Fq, Fr, G> element<C, Fq, Fr, G>::process_strauss_msm_rounds(const std::vector<element>& points,
                                                                        const std::vector<Fr>& scalars,
                                                                        const size_t max_num_bits)
{
    // Sanity checks
    BB_ASSERT_GT(points.size(), 0ULL, "process_strauss_msm: points cannot be empty");
    BB_ASSERT_EQ(points.size(), scalars.size(), "process_strauss_msm: points and scalars size mismatch");
 
    // Check that all scalars are in range
    for (const auto& scalar : scalars) {
        const size_t num_scalar_bits = static_cast<size_t>(uint512_t(scalar.get_value()).get_msb()) + 1ULL;
        BB_ASSERT_LTE(num_scalar_bits, max_num_bits, "process_strauss_msm: scalar out of range");
    }
 
    // Constant parameters
    const size_t num_rounds = max_num_bits;
    const size_t msm_size = scalars.size();
 
    // Compute ROM lookup table for points. Example if we have 3 points G1, G2, G3:
    // ┌───────┬─────────────────┐
    // │ Index │ Point           │
    // ├───────┼─────────────────┤
    // │   0   │  G1 + G2 + G3   │
    // │   1   │  G1 + G2 - G3   │
    // │   2   │  G1 - G2 + G3   │
    // │   3   │  G1 - G2 - G3   │
    // │   4   │ -G1 + G2 + G3   │
    // │   5   │ -G1 + G2 - G3   │
    // │   6   │ -G1 - G2 + G3   │
    // │   7   │ -G1 - G2 - G3   │
    // └───────┴─────────────────┘
    batch_lookup_table point_table(points);
 
    // Compute NAF representations of scalars
    std::vector<std::vector<bool_ct>> naf_entries;
    for (size_t i = 0; i < msm_size; ++i) {
        naf_entries.emplace_back(compute_naf(scalars[i], num_rounds));
    }
 
    // We choose a deterministic offset generator based on the number of rounds.
    // We compute both the initial and final offset generators: G_offset, 2ⁿ⁻¹ * G_offset.
    const auto [offset_generator_start, offset_generator_end] = compute_offset_generators(num_rounds);
 
    // Initialize accumulator with offset generator + first NAF column.
    element accumulator =
        element::chain_add_end(element::chain_add(offset_generator_start, point_table.get_chain_initial_entry()));
 
    // Process 4 NAF entries per iteration (for the remaining (num_rounds - 1) rounds)
    constexpr size_t num_rounds_per_iteration = 4;
    const size_t num_iterations = numeric::ceil_div((num_rounds - 1), num_rounds_per_iteration);
    const size_t num_rounds_per_final_iteration = (num_rounds - 1) - ((num_iterations - 1) * num_rounds_per_iteration);
 
    for (size_t i = 0; i < num_iterations; ++i) {
        std::vector<element::chain_add_accumulator> to_add;
        const size_t inner_num_rounds =
            (i != num_iterations - 1) ? num_rounds_per_iteration : num_rounds_per_final_iteration;
        for (size_t j = 0; j < inner_num_rounds; ++j) {
            // Gather the NAF columns for this iteration
            std::vector<bool_ct> nafs(msm_size);
            for (size_t k = 0; k < msm_size; ++k) {
                nafs[k] = (naf_entries[k][(i * num_rounds_per_iteration) + j + 1]);
            }
            to_add.emplace_back(point_table.get_chain_add_accumulator(nafs));
        }
 
        // Once we have looked-up all points from the four NAF columns, we update the accumulator as:
        // accumulator = 2.(2.(2.(2.accumulator + to_add[0]) + to_add[1]) + to_add[2]) + to_add[3]
        //             = 2⁴.accumulator + 2³.to_add[0] + 2².to_add[1] + 2¹.to_add[2] + to_add[3]
        accumulator = accumulator.multiple_montgomery_ladder(to_add);
    }
 
    // Subtract the skew factors (if any)
    for (size_t i = 0; i < msm_size; ++i) {
        element skew = accumulator - points[i];
        accumulator = accumulator.conditional_select(skew, naf_entries[i][num_rounds]);
    }
 
    // Subtract the scaled offset generator
    accumulator = accumulator - offset_generator_end;
 
    return accumulator;
}
 
template <typename C, class Fq, class Fr, class G>
element<C, Fq, Fr, G> element<C, Fq, Fr, G>::batch_mul(const std::vector<element>& _points,
                                                       const std::vector<Fr>& _scalars,
                                                       const size_t max_num_bits,
                                                       const bool with_edgecases,
                                                       const Fr& masking_scalar)
{
    // Sanity check input sizes
    BB_ASSERT_GT(_points.size(), 0ULL, "biggroup batch_mul: no points provided for batch multiplication");
    BB_ASSERT_EQ(_points.size(), _scalars.size(), "biggroup batch_mul: points and scalars size mismatch");
 
    // Replace (∞, scalar) pairs by the pair (G, 0).
    auto [points, scalars] = handle_points_at_infinity(_points, _scalars);
 
    BB_ASSERT_LTE(points.size(), _points.size());
    BB_ASSERT_EQ(points.size(),
                 scalars.size(),
                 "biggroup batch_mul: points and scalars size mismatch after handling points at infinity");
 
    // If batch_mul actually performs batch multiplication on the points and scalars, subprocedures can do
    // operations like addition or subtraction of points, which can trigger OriginTag security mechanisms
    // even though the final result satisfies the security logic. For example
    // result = submitted_in_round_0 * challenge_from_round_0 + submitted_in_round_1 * challenge_in_round_1
    // will trigger it, because the addition of submitted_in_round_0 to submitted_in_round_1 is dangerous by itself.
    // To avoid this, we remove the tags, merge them separately and set the result appropriately
    OriginTag tag{};
    const auto empty_tag = OriginTag();
    for (size_t i = 0; i < _points.size(); i++) {
        tag = OriginTag(tag, OriginTag(_points[i].get_origin_tag(), _scalars[i].get_origin_tag()));
    }
    for (size_t i = 0; i < scalars.size(); i++) {
        points[i].set_origin_tag(empty_tag);
        scalars[i].set_origin_tag(empty_tag);
    }
 
    // Accumulate constant-constant pairs out of circuit
    bool has_constant_terms = false;
    typename G::element constant_accumulator = G::element::infinity();
    std::vector<element> new_points;
    std::vector<Fr> new_scalars;
    for (size_t i = 0; i < points.size(); ++i) {
        if (points[i].is_constant() && scalars[i].is_constant()) {
            const auto& point_value = typename G::element(points[i].get_value());
            const auto& scalar_value = typename G::Fr(scalars[i].get_value());
            constant_accumulator += (point_value * scalar_value);
            has_constant_terms = true;
        } else {
            new_points.emplace_back(points[i]);
            new_scalars.emplace_back(scalars[i]);
        }
    }
    points = new_points;
    scalars = new_scalars;
 
    // If with_edgecases is false, masking_scalar must be constant and equal to 1 (as it is unused).
    if (!with_edgecases) {
        BB_ASSERT_EQ(
            masking_scalar.is_constant() && masking_scalar.get_value() == 1,
            true,
            "biggroup batch_mul: masking_scalar must be constant (and equal to 1) when with_edgecases is false");
    }
 
    if (with_edgecases) {
        // If points are linearly dependent, we randomise them using a masking scalar.
        // We do this to ensure that the x-coordinates of the points are all distinct. This is required
        // while creating the ROM lookup table with the points.
        std::tie(points, scalars) = mask_points(points, scalars, masking_scalar);
    }
 
    BB_ASSERT_EQ(
        points.size(), scalars.size(), "biggroup batch_mul: points and scalars size mismatch after handling edgecases");
 
    // Separate out zero scalars and corresponding points (because NAF(0) = NAF(modulus) which is 254 bits long)
    // Also add the last point and scalar to big_points and big_scalars (because its a 254-bit scalar)
    // We do this only if max_num_bits != 0 (i.e. we are not forced to use 254 bits anyway)
    const size_t original_size = scalars.size();
    std::vector<Fr> big_scalars;
    std::vector<element> big_points;
    std::vector<Fr> small_scalars;
    std::vector<element> small_points;
    for (size_t i = 0; i < original_size; ++i) {
        if (max_num_bits == 0) {
            big_points.emplace_back(points[i]);
            big_scalars.emplace_back(scalars[i]);
        } else {
            const bool is_last_scalar_big = ((i == original_size - 1) && with_edgecases);
            if (scalars[i].get_value() == 0 || is_last_scalar_big) {
                big_points.emplace_back(points[i]);
                big_scalars.emplace_back(scalars[i]);
            } else {
                small_points.emplace_back(points[i]);
                small_scalars.emplace_back(scalars[i]);
            }
        }
    }
 
    BB_ASSERT_EQ(original_size,
                 small_points.size() + big_points.size(),
                 "biggroup batch_mul: points size mismatch after separating big scalars");
    BB_ASSERT_EQ(big_points.size(),
                 big_scalars.size(),
                 "biggroup batch_mul: big points and scalars size mismatch after separating big scalars");
    BB_ASSERT_EQ(small_points.size(),
                 small_scalars.size(),
                 "biggroup batch_mul: small points and scalars size mismatch after separating big scalars");
 
    const size_t max_num_bits_in_field = Fr::modulus.get_msb() + 1;
 
    element accumulator;
    bool accumulator_initialized = false;
 
    // Check if we'll need to process any witness points
 
    // Initialize accumulator with constant terms if they exist, OR if there are no remaining points
    // (to handle the case where all points were filtered out by handle_points_at_infinity)
    const bool has_no_points = big_points.empty() && small_points.empty();
    if (has_constant_terms || has_no_points) {
        accumulator = element(constant_accumulator);
        accumulator_initialized = true;
    }
 
    if (!big_points.empty()) {
        // Process big scalars separately
        element big_result = element::process_strauss_msm_rounds(big_points, big_scalars, max_num_bits_in_field);
        accumulator = accumulator_initialized ? accumulator + big_result : big_result;
        accumulator_initialized = true;
    }
 
    if (!small_points.empty()) {
        // Process small scalars
        const size_t effective_max_num_bits = (max_num_bits == 0) ? max_num_bits_in_field : max_num_bits;
        element small_result = element::process_strauss_msm_rounds(small_points, small_scalars, effective_max_num_bits);
        accumulator = accumulator_initialized ? accumulator + small_result : small_result;
        accumulator_initialized = true;
    }
 
    accumulator.set_origin_tag(tag);
    return accumulator;
}
template <typename C, class Fq, class Fr, class G>
element<C, Fq, Fr, G> element<C, Fq, Fr, G>::operator*(const Fr& scalar) const
{
    // Use `scalar_mul` method without specifying the length of `scalar`.
    return scalar_mul(scalar);
}
 
template <typename C, class Fq, class Fr, class G>
element<C, Fq, Fr, G> element<C, Fq, Fr, G>::scalar_mul(const Fr& scalar, const size_t max_num_bits) const
{
    return element::batch_mul({ *this }, { scalar }, max_num_bits, /*with_edgecases=*/false);
}
} // namespace bb::stdlib::element_default