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/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)
    : x(x_in)
    , y(y_in)
    , _is_infinity(x.get_context() ? x.get_context() : y.get_context(), false)
{}
 
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)
    : x(x_in)
    , y(y_in)
    , _is_infinity(is_infinity)
{}
 
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^2 = x^3 + b.
    // If (x_1, y_1), (x_2, y_2) have x_1 == x_2, and the generic formula for lambda has a division by 0.
    // Then y_1 == y_2 (i.e. we are doubling) or y_2 == y_1 (the sum is infinity).
    // The cases have a special addition formula. 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 `lambda`. If we add, `lambda = (y2 - y1)/(x2 - x1)`, else `lambda = 3x1*x1/2y1
    const Fq add_lambda_numerator = other.y - y;
    const Fq xx = x * x;
    const Fq dbl_lambda_numerator = xx + xx + xx;
    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 inputs are points at infinity, we set lambda_denominator to be 1. This ensures we never trigger a
    // divide by zero error.
    // Note: if either inputs are points at infinity we will not use the result of this computation.
    Fq safe_edgecase_denominator = Fq(1);
    lambda_denominator =
        Fq::conditional_assign(has_infinity_input || infinity_predicate, safe_edgecase_denominator, lambda_denominator);
    const Fq lambda = Fq::div_without_denominator_check({ lambda_numerator }, lambda_denominator);
 
    const Fq x3 = lambda.sqradd({ -other.x, -x });
    const Fq y3 = lambda.madd(x - x3, { -y });
 
    element result(x3, y3);
    // 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);
 
    // is result point at infinity?
    // yes = infinity_predicate && !lhs_infinity && !rhs_infinity
    // yes = lhs_infinity && rhs_infinity
    // n.b. can likely optimize this
    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
{
 
    // if x_coordinates match, lambda triggers a divide by zero error.
    // Adding in `x_coordinates_match` ensures that lambda will always be well-formed
    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 `lambda`. If we add, `lambda = (y2 - y1)/(x2 - x1)`, else `lambda = 3x1*x1/2y1
    const Fq add_lambda_numerator = -other.y - y;
    const Fq xx = x * x;
    const Fq dbl_lambda_numerator = xx + xx + xx;
    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 inputs are points at infinity, we set lambda_denominator to be 1. This ensures we never trigger
    // a divide by zero error. (if either inputs are points at infinity we will not use the result of this
    // computation)
    Fq safe_edgecase_denominator = Fq(1);
    lambda_denominator =
        Fq::conditional_assign(has_infinity_input || infinity_predicate, safe_edgecase_denominator, lambda_denominator);
    const Fq lambda = Fq::div_without_denominator_check({ lambda_numerator }, lambda_denominator);
 
    const Fq x3 = lambda.sqradd({ -other.x, -x });
    const Fq y3 = lambda.madd(x - x3, { -y });
 
    element result(x3, y3);
    // 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);
 
    // is result point at infinity?
    // yes = infinity_predicate && !lhs_infinity && !rhs_infinity
    // yes = lhs_infinity && rhs_infinity
    // n.b. can likely optimize this
    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);
}
 
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);
}
 
// TODO(https://github.com/AztecProtocol/barretenberg/issues/657): This function is untested
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
{
 
    // TODO(https://github.com/AztecProtocol/barretenberg/issues/971): This will fail when the two elements are
    // the same even in the case of a valid circuit
    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), element(x_4, y_4) };
}
 
template <typename C, class Fq, class Fr, class G> element<C, Fq, Fr, G> element<C, Fq, Fr, G>::dbl() const
{
 
    Fq two_x = x + x;
    if constexpr (G::has_a) {
        Fq a(get_context(), uint256_t(G::curve_a));
        Fq neg_lambda = Fq::msub_div({ x }, { (two_x + x) }, (y + y), { a }, /*enable_divisor_nz_check*/ false);
        Fq x_3 = neg_lambda.sqradd({ -(two_x) });
        Fq y_3 = neg_lambda.madd(x_3 - x, { -y });
        // TODO(suyash): do we handle the point at infinity case here?
        return element(x_3, y_3);
    }
    // TODO(): handle y = 0 case.
    Fq neg_lambda = Fq::msub_div({ x }, { (two_x + x) }, (y + y), {}, /*enable_divisor_nz_check*/ false);
    Fq x_3 = neg_lambda.sqradd({ -(two_x) });
    Fq y_3 = neg_lambda.madd(x_3 - x, { -y });
    element result = element(x_3, y_3);
    result.set_point_at_infinity(is_point_at_infinity());
    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;
 
    p1.x.assert_is_not_equal(p2.x);
    const Fq lambda = Fq::div_without_denominator_check({ p2.y, -p1.y }, (p2.x - p1.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)
{
    // use `chain_add_start` to start an addition chain (i.e. if acc has a y-coordinate)
    if (acc.is_element) {
        return chain_add_start(p1, element(acc.x3_prev, acc.y3_prev));
    }
    // validate we can use incomplete addition formulae
    p1.x.assert_is_not_equal(acc.x3_prev);
 
    // lambda = (y2 - y1) / (x2 - x1)
    // but we don't have y2!
    // however, we do know that y2 = lambda_prev * (x1_prev - x2) - y1_prev
    // => lambda * (x2 - x1) = lambda_prev * (x1_prev - x2) - y1_prev - y1
    // => lambda * (x2 - x1) + lambda_prev * (x2 - x1_prev) + y1 + y1_pev = 0
    // => lambda = lambda_prev * (x1_prev - x2) - y1_prev - y1 / (x2 - x1)
    // => lambda = - (lambda_prev * (x2 - x1_prev) + y1_prev + y1) / (x2 - x1)
 
    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 x2 != p1.x is enforced
    const auto x3 = lambda.sqradd({ -x2, -p1.x });
 
    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 (acc.is_element) {
        return element(acc.x3_prev, acc.y3_prev);
    }
    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);
}
// #################################
// ### SCALAR MULTIPLICATION METHODS
// #################################
template <typename C, class Fq, class Fr, class G>
element<C, Fq, Fr, G> element<C, Fq, Fr, G>::montgomery_ladder(const element& other) const
{
    other.x.assert_is_not_equal(x);
    const Fq lambda_1 = Fq::div_without_denominator_check({ other.y - y }, (other.x - x));
 
    const Fq x_3 = lambda_1.sqradd({ -other.x, -x });
 
    const Fq minus_lambda_2 = lambda_1 + Fq::div_without_denominator_check({ y + y }, (x_3 - x));
 
    const Fq x_4 = minus_lambda_2.sqradd({ -x, -x_3 });
 
    const Fq y_4 = minus_lambda_2.madd(x_4 - x, { -y });
    return element(x_4, y_4);
}
 
template <typename C, class Fq, class Fr, class G>
element<C, Fq, Fr, G> element<C, Fq, Fr, G>::montgomery_ladder(const chain_add_accumulator& to_add)
{
    if (to_add.is_element) {
        throw_or_abort("An accumulator expected");
    }
    x.assert_is_not_equal(to_add.x3_prev);
 
    // lambda = (y2 - y1) / (x2 - x1)
    // but we don't have y2!
    // however, we do know that y2 = lambda_prev * (x1_prev - x2) - y1_prev
    // => lambda * (x2 - x1) = lambda_prev * (x1_prev - x2) - y1_prev - y1
    // => lambda * (x2 - x1) + lambda_prev * (x2 - x1_prev) + y1 + y1_pev = 0
    // => lambda = lambda_prev * (x1_prev - x2) - y1_prev - y1 / (x2 - x1)
    // => lambda = - (lambda_prev * (x2 - x1_prev) + y1_prev + y1) / (x2 - x1)
 
    auto& x2 = to_add.x3_prev;
    const auto lambda = Fq::msub_div({ to_add.lambda_prev },
                                     { (x2 - to_add.x1_prev) },
                                     (x2 - x),
                                     { to_add.y1_prev, y },
                                     /*enable_divisor_nz_check*/ false); // divisor is non-zero as x2 != x is enforced
    const auto x3 = lambda.sqradd({ -x2, -x });
 
    const Fq minus_lambda_2 = lambda + Fq::div_without_denominator_check({ y + y }, (x3 - x));
 
    const Fq x4 = minus_lambda_2.sqradd({ -x, -x3 });
 
    const Fq y4 = minus_lambda_2.madd(x4 - x, { -y });
    return element(x4, y4);
}
 
template <typename C, class Fq, class Fr, class G>
element<C, Fq, Fr, G> element<C, Fq, Fr, G>::quadruple_and_add(const std::vector<element>& to_add) const
{
    const Fq two_x = x + x;
    Fq x_1;
    Fq minus_lambda_dbl;
    if constexpr (G::has_a) {
        Fq a(get_context(), uint256_t(G::curve_a));
        minus_lambda_dbl = Fq::msub_div({ x }, { (two_x + x) }, (y + y), { a }, /*enable_divisor_nz_check*/ false);
        x_1 = minus_lambda_dbl.sqradd({ -(two_x) });
    } else {
        // TODO(): handle y = 0 case.
        minus_lambda_dbl = Fq::msub_div({ x }, { (two_x + x) }, (y + y), {}, /*enable_divisor_nz_check*/ false);
        x_1 = minus_lambda_dbl.sqradd({ -(two_x) });
    }
 
    BB_ASSERT_GT(to_add.size(), 0);
    to_add[0].x.assert_is_not_equal(x_1);
 
    const Fq x_minus_x_1 = x - x_1;
 
    const Fq lambda_1 =
        Fq::msub_div({ minus_lambda_dbl },
                     { x_minus_x_1 },
                     (x_1 - to_add[0].x),
                     { to_add[0].y, y },
                     /*enable_divisor_nz_check*/ false); // divisor is non-zero as x1 != to_add[0].x is enforced
 
    const Fq x_3 = lambda_1.sqradd({ -to_add[0].x, -x_1 });
 
    // TODO(): analyse if (x3 - x1) can be zero.
    const Fq half_minus_lambda_2_minus_lambda_1 =
        Fq::msub_div({ minus_lambda_dbl },
                     { x_minus_x_1 },
                     (x_3 - x_1),
                     { y },
                     /*enable_divisor_nz_check*/ false); // divisor is non-zero as x3 != x1 is enforced
 
    const Fq minus_lambda_2_minus_lambda_1 = half_minus_lambda_2_minus_lambda_1 + half_minus_lambda_2_minus_lambda_1;
    const Fq minus_lambda_2 = minus_lambda_2_minus_lambda_1 + lambda_1;
 
    const Fq x_4 = minus_lambda_2.sqradd({ -x_1, -x_3 });
 
    const Fq x_4_sub_x_1 = x_4 - x_1;
 
    if (to_add.size() == 1) {
        const Fq y_4 = Fq::dual_madd(minus_lambda_2, x_4_sub_x_1, minus_lambda_dbl, x_minus_x_1, { y });
        return element(x_4, y_4);
    }
    to_add[1].x.assert_is_not_equal(to_add[0].x);
 
    // TODO(): analyse if (x4 - x1) can be zero.
    Fq minus_lambda_3 = Fq::msub_div({ minus_lambda_dbl, minus_lambda_2 },
                                     { x_minus_x_1, x_4_sub_x_1 },
                                     (x_4 - to_add[1].x),
                                     { y, -(to_add[1].y) },
                                     /*enable_divisor_nz_check*/ false);
 
    // X5 = L3.L3 - X4 - XB
    const Fq x_5 = minus_lambda_3.sqradd({ -x_4, -to_add[1].x });
 
    if (to_add.size() == 2) {
        // Y5 = L3.(XB - X5) - YB
        const Fq y_5 = minus_lambda_3.madd(x_5 - to_add[1].x, { -to_add[1].y });
        return element(x_5, y_5);
    }
 
    Fq x_prev = x_5;
    Fq minus_lambda_prev = minus_lambda_3;
 
    for (size_t i = 2; i < to_add.size(); ++i) {
 
        to_add[i].x.assert_is_not_equal(to_add[i - 1].x);
        // Lambda = Yprev - Yadd[i] / Xprev - Xadd[i]
        //        = -Lprev.(Xprev - Xadd[i-1]) - Yadd[i - 1] - Yadd[i] / Xprev - Xadd[i]
        const Fq minus_lambda =
            Fq::msub_div({ minus_lambda_prev },
                         { to_add[i - 1].x - x_prev },
                         (to_add[i].x - x_prev),
                         { to_add[i - 1].y, to_add[i].y },
                         /*enable_divisor_nz_check*/ false); // divisor is non-zero as x_prev != to_add[i].x is enforced
 
        // X = Lambda * Lambda - Xprev - Xadd[i]
        const Fq x_out = minus_lambda.sqradd({ -x_prev, -to_add[i].x });
 
        x_prev = x_out;
        minus_lambda_prev = minus_lambda;
    }
    const Fq y_out = minus_lambda_prev.madd(x_prev - to_add[to_add.size() - 1].x, { -to_add[to_add.size() - 1].y });
 
    return element(x_prev, y_out);
}
 
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;
    };
 
    Fq previous_x = x;
    composite_y previous_y{ std::vector<Fq>(), std::vector<Fq>(), std::vector<Fq>(), false };
    for (size_t i = 0; i < add.size(); ++i) {
        previous_x.assert_is_not_equal(add[i].x3_prev);
 
        // composite_y add_y;
        bool negate_add_y = (i > 0) && !previous_y.is_negative;
        std::vector<Fq> lambda1_left;
        std::vector<Fq> lambda1_right;
        std::vector<Fq> lambda1_add;
 
        if (i == 0) {
            lambda1_add.emplace_back(-y);
        } else {
            lambda1_left = previous_y.mul_left;
            lambda1_right = previous_y.mul_right;
            lambda1_add = previous_y.add;
        }
 
        if (!add[i].is_element) {
            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 if (i > 0) {
            lambda1_add.emplace_back(negate_add_y ? -add[i].y3_prev : add[i].y3_prev);
        }
        // if previous_y is negated then add stays positive
        // if previous_y is positive then add stays negated
        // | add.y is negated | previous_y is negated | output of msub_div is -lambda |
        // | --- | --- | --- |
        // | no  | yes | yes |
        // | yes | no  | no  |
 
        Fq lambda1;
        if (!add[i].is_element || i > 0) {
            bool flip_lambda1_denominator = !negate_add_y;
            Fq denominator = flip_lambda1_denominator ? previous_x - add[i].x3_prev : add[i].x3_prev - previous_x;
            lambda1 = Fq::msub_div(
                lambda1_left,
                lambda1_right,
                denominator,
                lambda1_add,
                /*enable_divisor_nz_check*/ false); // divisor is non-zero as previous_x != add[i].x3_prev is enforced
        } else {
            lambda1 = Fq::div_without_denominator_check({ add[i].y3_prev - y }, (add[i].x3_prev - x));
        }
 
        Fq x_3 = lambda1.madd(lambda1, { -add[i].x3_prev, -previous_x });
 
        // We can avoid computing y_4, instead substituting the expression `minus_lambda_2 * (x_4 - x) - y`
        // where needed. This is cheaper, because we can evaluate two field multiplications (or a field
        // multiplication + a field division) with only one non-native field reduction. E.g. evaluating (a * b)
        // + (c * d) = e mod p only requires 1 quotient and remainder, which is the major cost of a non-native
        // field multiplication
        Fq lambda2;
        if (i == 0) {
            lambda2 = Fq::div_without_denominator_check({ y + y }, (previous_x - x_3)) - lambda1;
        } else {
            Fq l2_denominator = previous_y.is_negative ? previous_x - x_3 : x_3 - previous_x;
            // TODO(): analyse if l2_denominator can be zero.
            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;
        }
 
        Fq x_4 = lambda2.sqradd({ -x_3, -previous_x });
        composite_y y_4;
        if (i == 0) {
            // We want to make sure that at the final iteration, `y_previous.is_negative = false`
            // Each iteration flips the sign of y_previous.is_negative.
            // i.e. whether we store y_4 or -y_4 depends on the number of points we have
            bool num_points_even = ((add.size() & 0x01UL) == 0);
            y_4.add.emplace_back(num_points_even ? y : -y);
            y_4.mul_left.emplace_back(lambda2);
            y_4.mul_right.emplace_back(num_points_even ? x_4 - previous_x : previous_x - x_4);
            y_4.is_negative = num_points_even;
        } else {
            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 in 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;
 
    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);
}
 
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>::batch_mul(const std::vector<element>& _points,
                                                       const std::vector<Fr>& _scalars,
                                                       const size_t max_num_bits,
                                                       const bool with_edgecases)
{
    auto [points, scalars] = handle_points_at_infinity(_points, _scalars);
    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++) {
        // 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
        points[i].set_origin_tag(empty_tag);
        scalars[i].set_origin_tag(empty_tag);
    }
 
    // Perform goblinized batched mul if available; supported only for BN254
    if (with_edgecases) {
        std::tie(points, scalars) = mask_points(points, scalars);
    }
    const size_t num_points = points.size();
    BB_ASSERT_EQ(scalars.size(), num_points);
 
    batch_lookup_table point_table(points);
    const size_t num_rounds = (max_num_bits == 0) ? Fr::modulus.get_msb() + 1 : max_num_bits;
 
    std::vector<std::vector<bool_ct>> naf_entries;
    for (size_t i = 0; i < num_points; ++i) {
        naf_entries.emplace_back(compute_naf(scalars[i], max_num_bits));
    }
    const auto offset_generators = compute_offset_generators(num_rounds);
    element accumulator =
        element::chain_add_end(element::chain_add(offset_generators.first, point_table.get_chain_initial_entry()));
 
    constexpr size_t num_rounds_per_iteration = 4;
    size_t num_iterations = num_rounds / num_rounds_per_iteration;
    num_iterations += ((num_iterations * num_rounds_per_iteration) == num_rounds) ? 0 : 1;
    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<bool_ct> nafs(num_points);
        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) {
            for (size_t k = 0; k < num_points; ++k) {
                nafs[k] = (naf_entries[k][i * num_rounds_per_iteration + j + 1]);
            }
            to_add.emplace_back(point_table.get_chain_add_accumulator(nafs));
        }
        accumulator = accumulator.multiple_montgomery_ladder(to_add);
    }
    for (size_t i = 0; i < num_points; ++i) {
        element skew = accumulator - points[i];
        Fq out_x = accumulator.x.conditional_select(skew.x, naf_entries[i][num_rounds]);
        Fq out_y = accumulator.y.conditional_select(skew.y, naf_entries[i][num_rounds]);
        accumulator = element(out_x, out_y);
    }
    accumulator = accumulator - offset_generators.second;
 
    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
{
    BB_ASSERT_EQ(max_num_bits % 2, 0U);
    OriginTag tag{};
    tag = OriginTag(tag, OriginTag(this->get_origin_tag(), scalar.get_origin_tag()));
 
    bool_ct is_point_at_infinity = this->is_point_at_infinity();
 
    const size_t num_rounds = (max_num_bits == 0) ? Fr::modulus.get_msb() + 1 : max_num_bits;
 
    element result;
    if (max_num_bits != 0) {
        // The case of short scalars
        result = element::bn254_endo_batch_mul({}, {}, { *this }, { scalar }, num_rounds);
    } else {
        // The case of arbitrary length scalars
        result = element::bn254_endo_batch_mul({ *this }, { scalar }, {}, {}, num_rounds);
    };
 
    // Handle point at infinity
    result.x = Fq::conditional_assign(is_point_at_infinity, x, result.x);
    result.y = Fq::conditional_assign(is_point_at_infinity, y, result.y);
 
    result.set_point_at_infinity(is_point_at_infinity);
 
    // Propagate the origin tag
    result.set_origin_tag(tag);
 
    return result;
}
} // namespace bb::stdlib::element_default