biggroup_impl.hpp

Go to the documentation of this file.

// === AUDIT STATUS ===
// internal:    { status: Complete, auditors: [Suyash], commit: 553c5eb82901955c638b943065acd3e47fc918c0}
// external_1:  { status: not started, auditors: [], commit: }
// external_2:  { status: not started, auditors: [], commit: }
// =====================
 
#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 Fq& x_in, const Fq& y_in, const bool assert_on_curve)
    : _x(x_in)
    , _y(y_in)
{
    // Detect infinity: point is at infinity iff both coordinates are zero.
    // We sum all 8 binary basis limbs (4 from x, 4 from y) and check if the sum is zero.
    // This works because: (1) after reduction, each limb is non-negative and range-constrained,
    // so sum=0 iff all limbs=0; (2) max sum is 8 * 2^68 ≈ 2^71 << n ≈ 2^254, so no native wraparound.
    field_ct limb_sum = 0;
    for (size_t i = 0; i < Fq::NUM_LIMBS; ++i) {
        limb_sum += _x.get_limb(i).element;
    }
    for (size_t i = 0; i < Fq::NUM_LIMBS; ++i) {
        limb_sum += _y.get_limb(i).element;
    }
    _is_infinity = limb_sum.is_zero();
 
    // Validate on-curve if requested
    if (assert_on_curve) {
        [[maybe_unused]] Fq _ = 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 stdlib::bool_t<C>& is_infinity,
                               const bool assert_on_curve)
    : _x(x_in)
    , _y(y_in)
    , _is_infinity(is_infinity.normalize())
{
    BB_ASSERT(!(_x.is_constant() && _y.is_constant() && !_is_infinity.is_constant()),
              "biggroup: constant coordinates with non-constant infinity flag");
    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>::add_internal(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.
    // Exception: Points at infinity use the canonical (0, 0) representation, so skip the check for them.
    if (!lhs_infinity.get_value()) {
        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.");
    }
    if (!rhs_infinity.get_value()) {
        const typename G::Fq other_y_value = uint256_t(other._y.get_value());
        BB_ASSERT_EQ((other_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);
 
    // Compute λ = numerator / denominator
    // We enforce that denominator is not zero to protect against soundness issues.
    const Fq lambda = lambda_numerator / lambda_denominator;
 
    // Compute resulting point coordinates: x₃ = λ² - x₁ - x₂, y₃ = λ(x₁ - x₃) - y₁
    Fq x3 = lambda.sqradd({ -other._x, -_x });
    Fq y3 = lambda.madd(_x - x3, { -_y });
 
    // if lhs infinity, return rhs
    x3 = Fq::conditional_assign(lhs_infinity, other._x, x3);
    y3 = Fq::conditional_assign(lhs_infinity, other._y, y3);
    // if rhs infinity, return lhs
    x3 = Fq::conditional_assign(rhs_infinity, _x, x3);
    y3 = Fq::conditional_assign(rhs_infinity, _y, y3);
 
    // 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);
    mark_witness_as_used(field_t<C>(result_is_infinity));
 
    element result(x3, y3, /*is_infinity=*/result_is_infinity, /*assert_on_curve=*/false);
    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>::operator+(const element& other) const
{
    element result = add_internal(other);
    return result.get_standard_form();
}
 
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(get_context(), 0);
    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>::subtract_internal(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.
    // Defense in depth: also guard when doubling with y = 0 (denominator 2y would be zero).
    const bool_ct safe_denominator_needed = has_infinity_input || infinity_predicate;
    lambda_denominator = Fq::conditional_assign(safe_denominator_needed, Fq(1), lambda_denominator);
 
    // Now compute lambda = numerator / denominator
    // We enforce that the denominator is not zero (in division operator), so this division is safe.
    const Fq lambda = lambda_numerator / lambda_denominator;
 
    // Compute resulting point coordinates: x₃ = λ² - x₁ - x₂, y₃ = λ(x₁ - x₃) - y₁
    Fq x3 = lambda.sqradd({ -other._x, -_x });
    Fq y3 = lambda.madd(_x - x3, { -_y });
 
    // if lhs infinity, return -rhs (negated rhs point)
    x3 = Fq::conditional_assign(lhs_infinity, other._x, x3);
    y3 = Fq::conditional_assign(lhs_infinity, -other._y, y3);
    // if rhs infinity, return lhs
    x3 = Fq::conditional_assign(rhs_infinity, _x, x3);
    y3 = Fq::conditional_assign(rhs_infinity, _y, y3);
 
    // 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);
    mark_witness_as_used(field_t<C>(result_is_infinity));
 
    element result(x3, y3, /*is_infinity=*/result_is_infinity, /*assert_on_curve=*/false);
    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>::operator-(const element& other) const
{
    element result = subtract_internal(other);
    return result.get_standard_form();
}
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 });
    // Use 4-arg constructor with is_infinity=false (checked operations assume valid, non-infinity points).
    // This avoids the expensive bigfield equality checks in the 2-arg constructor's infinity auto-detection.
    return element(x3, y3, bool_ct(x3.get_context(), false), /*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 });
 
    // Use 4-arg constructor with is_infinity=false (checked operations assume valid, non-infinity points).
    return element(x_3, y_3, bool_ct(x_3.get_context(), false), /*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 });
 
    // Use 4-arg constructor with is_infinity=false (checked operations assume valid, non-infinity points).
    bool_ct not_infinity(x_3.get_context(), false);
    return { element(x_3, y_3, not_infinity, /*assert_on_curve=*/false),
             element(x_4, y_4, not_infinity, /*assert_on_curve=*/false) };
}
 
template <typename C, class Fq, class Fr, class G> element<C, Fq, Fr, G> element<C, Fq, Fr, G>::dbl_internal() const
{
    const bool_ct is_infinity = is_point_at_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.
    // Exception: Points at infinity use the canonical (0, 0) representation, so skip the check for them.
    if (!is_infinity.get_value()) {
        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.");
    }
 
    // If the input is a point at infinity, use a safe denominator (1) to prevent division by zero.
    // The result will be infinity anyway, so the computed coordinates don't matter.
    const Fq two_y = _y + _y;
    Fq denominator = Fq::conditional_assign(is_infinity, Fq(get_context(), 1), two_y);
 
    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)
        // We enforce the denominator is not zero to protect against soundness issues.
        Fq neg_lambda = Fq::msub_div({ _x }, { (two_x + _x) }, denominator, { a }, /*enable_divisor_nz_check*/ true);
 
        // 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 });
 
        mark_witness_as_used(field_t<C>(is_infinity));
        return element(x_3, y_3, /*is_infinity=*/is_infinity, /*assert_on_curve=*/false);
    }
 
    // Curve equation when a = 0: y² = x³ + b
    // Compute neg_lambda = -λ = -3x² / (2y)
    // msub_div computes: -(Σᵢ aᵢ·bᵢ) / d = -(x·(3x)) / (2y) = -3x² / (2y)
    // We enforce the denominator is not zero to protect against soundness issues.
    Fq neg_lambda = Fq::msub_div({ _x }, { (two_x + _x) }, denominator, {}, /*enable_divisor_nz_check*/ true);
 
    // 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 });
 
    mark_witness_as_used(field_t<C>(is_infinity));
    return element(x_3, y_3, /*is_infinity=*/is_infinity, /*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
{
    element result = dbl_internal();
    return result.get_standard_form();
}
 
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) {
        // Use 4-arg constructor with is_infinity=false (chain operations assume valid, non-infinity points).
        return chain_add_start(
            p1,
            element(acc.x3_prev, acc.y3_prev, bool_ct(acc.x3_prev.get_context(), false), /*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) {
        // Use 4-arg constructor with is_infinity=false (chain operations assume valid, non-infinity points).
        return element(acc.x3_prev, acc.y3_prev, bool_ct(acc.x3_prev.get_context(), false), /*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 });
    // Use 4-arg constructor with is_infinity=false (chain operations assume valid, non-infinity points).
    return element(x3, y3, bool_ct(x3.get_context(), false), /*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, "biggroup::multiple_montgomery_ladder: x-coordinates must be distinct.");
 
    // 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, "biggroup::multiple_montgomery_ladder: x-coordinates must be distinct.");
    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,
                                       "biggroup::multiple_montgomery_ladder: x-coordinates must be distinct.");
 
        // 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, "biggroup::multiple_montgomery_ladder: x-coordinates must be distinct.");
        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)
        y_4.mul_left.insert(y_4.mul_left.end(), previous_y.mul_left.begin(), previous_y.mul_left.end());
        y_4.mul_right.insert(y_4.mul_right.end(), previous_y.mul_right.begin(), previous_y.mul_right.end());
        y_4.add.insert(y_4.add.end(), previous_y.add.begin(), previous_y.add.end());
 
        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);
    // Use 4-arg constructor with is_infinity=false (montgomery ladder assumes valid, non-infinity points).
    return element(x_out, y_out, bool_ct(x_out.get_context(), false), /*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>(element(offset_generator.x, offset_generator.y),
                                            element(offset_generator_end.x, offset_generator_end.y));
}
 
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.subtract_internal(points[i]);
        accumulator = accumulator.conditional_select(skew, naf_entries[i][num_rounds]);
    }
 
    // Subtract the scaled offset generator
    accumulator = accumulator.subtract_internal(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_internal(const std::vector<element>& _points,
                                                                const std::vector<Fr>& _scalars,
                                                                const size_t max_num_bits,
                                                                const bool with_edgecases)
{
    // 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");
 
    // Get builder context from input points (needed for creating infinity results)
    C* builder = _points[0].get_context();
 
    // 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 = OriginTag::constant();  // Initialize as CONSTANT so merging with input tags works correctly
    auto empty_tag = OriginTag::constant(); // Disable origin checking during intermediate operations
    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 && !points.empty()) {
        // If points are linearly dependent, we randomise them using a free-witness offset generator.
        // 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.
        auto [masked_points, masked_scalars, _offset_generator] = mask_points(points, scalars);
        points = std::move(masked_points);
        scalars = std::move(masked_scalars);
    }
 
    BB_ASSERT_EQ(
        points.size(), scalars.size(), "biggroup batch_mul: points and scalars size mismatch after handling edgecases");
 
    // Partition scalars into big and small paths using compile-time information.
    // Since `with_edgecases` appends a 254-bit randomization scalar at the end (via
    // `mask_points`), the last scalar must go to big. Everything else routes by
    // `max_num_bits`: `0` forces the full-width path, anything else uses the small path.
    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) {
        const bool is_last_scalar_big = (i == original_size - 1) && with_edgecases;
        if (max_num_bits == 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) {
        // Convert from projective to affine to get correct (x, y) coordinates
        typename G::affine_element constant_accumulator_affine(constant_accumulator);
        if (constant_accumulator_affine.is_point_at_infinity()) {
            // For infinity, create element with canonical (0, 0) coordinates
            // Using constant zero Fq to create a constant infinity element
            Fq zero_fq = Fq(builder, 0);
            accumulator = element(zero_fq, zero_fq);
        } else {
            accumulator = element(constant_accumulator_affine.x, constant_accumulator_affine.y);
        }
        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.add_internal(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.add_internal(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>::batch_mul(const std::vector<element>& points,
                                                       const std::vector<Fr>& scalars,
                                                       const size_t max_num_bits,
                                                       const bool with_edgecases)
{
    element result = batch_mul_internal(points, scalars, max_num_bits, with_edgecases);
    return result.get_standard_form();
}
 
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