Optimize ecmul and ecrecovery implementation with curve endomorphism

lib/evmone_precompiles/bn254.cpp

@@ -12,6 +12,28 @@ namespace
 const ModArith<uint256> Fp{FieldPrime};
 const auto B = Fp.to_mont(3);
 const auto B3 = Fp.to_mont(3 * 3);
+
+struct Config
+{
+    // Linearly independent short vectors (𝑣₁=(𝑥₁, 𝑦₁), 𝑣₂=(x₂, 𝑦₂)) such that f(𝑣₁) = f(𝑣₂) = 0,
+    // where f : ℤ×ℤ → ℤₙ is defined as (𝑖,𝑗) → (𝑖+𝑗λ), where λ² + λ ≡ -1 mod n. n is bn245 curve
+    // order. Here λ = 0xb3c4d79d41a917585bfc41088d8daaa78b17ea66b99c90dd. DET is (𝑣₁, 𝑣₂) matrix
+    // determinant. For more details see https://www.iacr.org/archive/crypto2001/21390189.pdf
+    static constexpr auto X1 = 147946756881789319020627676272574806254_u512;
+    // Y1 should be negative, hence we calculate the determinant below adding operands instead of
+    // subtracting.
+    static constexpr auto Y1 = 147946756881789318990833708069417712965_u512;
+    static constexpr auto X2 = 147946756881789319000765030803803410728_u512;
+    static constexpr auto Y2 = 147946756881789319010696353538189108491_u512;
+    static constexpr auto DET =
+        43776485743678550444492811490514550177096728800832068687396408373151616991234_u256;
+};
+
+// For bn254 curve and β ∈ 𝔽ₚ endomorphism ϕ : E₂ → E₂ defined as (𝑥,𝑦) → (β𝑥,𝑦) calculates [λ](𝑥,𝑦)
+// with only one multiplication in 𝔽ₚ. BETA value in Montgomery form;
+inline constexpr auto BETA =
+    20006444479023397533370224967097343182639219473961804911780625968796493078869_u256;
+
 }  // namespace
 
 bool validate(const Point& pt) noexcept
@@ -49,9 +71,16 @@ Point mul(const Point& pt, const uint256& c) noexcept
     if (c == 0)
         return {};
 
-    const auto pr = ecc::mul(Fp, ecc::to_proj(Fp, pt), c, B3);
+    const auto [k1, k2] = ecc::decompose<Config>(c % Order);
 
-    return ecc::to_affine(Fp, field_inv, pr);
+    const ecc::ProjPoint<uint256> q = {Fp.mul(BETA, Fp.to_mont(pt.x)),
+        !k2.first ? Fp.to_mont(pt.y) : Fp.to_mont(FieldPrime - pt.y), Fp.to_mont(1)};
+
+    const auto r = shamir_multiply(Fp, B3, k1.second,
+        !k1.first ? ecc::to_proj(Fp, pt) : ecc::to_proj(Fp, {pt.x, FieldPrime - pt.y}), k2.second,
+        q);
+
+    return ecc::to_affine(Fp, field_inv, r);
 }
 
 uint256 field_inv(const ModArith<uint256>& m, const uint256& x) noexcept

lib/evmone_precompiles/bn254.hpp

@@ -13,6 +13,10 @@ using namespace intx;
 inline constexpr auto FieldPrime =
     0x30644e72e131a029b85045b68181585d97816a916871ca8d3c208c16d87cfd47_u256;
 
+/// The bn254 curve order (N)
+inline constexpr auto Order =
+    0x30644E72E131A029B85045B68181585D2833E84879B9709143E1F593F0000001_u256;
+
 using Point = ecc::Point<uint256>;
 
 /// Validates that point is from the bn254 curve group

lib/evmone_precompiles/ecc.hpp

@@ -202,5 +202,83 @@ ProjPoint<IntT> mul(const evmmax::ModArith<IntT>& s, const ProjPoint<IntT>& z, c
     return p;
 }
 
+// Computes uG + vQ using "Shamir's trick". https://eprint.iacr.org/2003/257.pdf (page 7)
+template <typename UIntT>
+inline ProjPoint<UIntT> shamir_multiply(const ModArith<UIntT>& m, const UIntT& b3, const UIntT& u,
+    const ProjPoint<UIntT>& g, const UIntT& v, const ProjPoint<UIntT>& q)
+{
+    ProjPoint<UIntT> r;
+    const auto h = add(m, g, q, b3);
+
+    const auto u_lz = clz(u);
+    const auto v_lz = clz(v);
+
+    auto lz = std::min(u_lz, v_lz);
+
+    if (lz == UIntT::num_bits)
+        return {};
+
+    if (u_lz < v_lz)
+        r = g;
+    else if (u_lz > v_lz)
+        r = q;
+    else
+        r = h;
+
+    auto mask = (UIntT{1} << (UIntT::num_bits - 1 - lz - 1));
+
+    while (mask != 0)
+    {
+        r = dbl(m, r, b3);
+        if (u & v & mask)
+            r = add(m, r, h, b3);
+        else if (u & mask)
+            r = add(m, r, g, b3);
+        else if (v & mask)
+            r = add(m, r, q, b3);
+
+        mask >>= 1;
+    }
+
+    return r;
+}
+
+// Decomposes scalar k into k₁ and k₂ such that k₁ + k₂λ ≡ k mod n
+// Returns ((is_negative, k1), (is_negative, k2))
+template <typename ConfigT, typename UIntT>
+inline std::pair<std::pair<bool, UIntT>, std::pair<bool, UIntT>> decompose(const UIntT& k) noexcept
+{
+    using DIntT = intx::uint<2 * UIntT::num_bits>;
+
+    const auto round_div = [](const DIntT& n) {
+        const auto [q, r] = udivrem(n, ConfigT::DET);
+
+        return (r <= (ConfigT::DET / 2)) ? q : (q + 1);
+    };
+
+    const auto z1 = round_div(ConfigT::Y2 * k);
+    const auto z2 = round_div(ConfigT::Y1 * k);
+
+    auto const z1x1_z2x2 = z1 * ConfigT::X1 + z2 * ConfigT::X2;
+
+    auto k1_is_neg = false;
+    auto k2_is_neg = false;
+
+    auto tk = k;
+    if (tk < z1x1_z2x2)
+        k1_is_neg = true;
+
+    const auto k1 = !k1_is_neg ? (tk - z1x1_z2x2) : z1x1_z2x2 - tk;
+
+    const DIntT z2y2 = z2 * ConfigT::Y2;
+    const DIntT z1y1 = z1 * ConfigT::Y1;
+
+    if (z1y1 < z2y2)
+        k2_is_neg = true;
+
+    const DIntT k2 = !k2_is_neg ? (z1y1 - z2y2) : z2y2 - z1y1;
+
+    return {{k1_is_neg, UIntT{k1}}, {k2_is_neg, UIntT{k2}}};
+}
 
 }  // namespace evmmax::ecc

lib/evmone_precompiles/secp256k1.cpp

@@ -3,6 +3,7 @@
 // SPDX-License-Identifier: Apache-2.0
 #include "secp256k1.hpp"
 #include <ethash/keccak.hpp>
+#include <iostream>
 
 namespace evmmax::secp256k1
 {
@@ -14,6 +15,28 @@ const auto B3 = Fp.to_mont(7 * 3);
 
 constexpr Point G{0x79be667ef9dcbbac55a06295ce870b07029bfcdb2dce28d959f2815b16f81798_u256,
     0x483ada7726a3c4655da4fbfc0e1108a8fd17b448a68554199c47d08ffb10d4b8_u256};
+
+struct Config
+{
+    // Linearly independent short vectors (𝑣₁=(𝑥₁, 𝑦₁), 𝑣₂=(x₂, 𝑦₂)) such that f(𝑣₁) = f(𝑣₂) = 0,
+    // where f : ℤ×ℤ → ℤₙ is defined as (𝑖,𝑗) → (𝑖+𝑗λ), where λ² + λ ≡ -1 mod n. n is secp256k1
+    // curve order. Here λ = 0x5363ad4cc05c30e0a5261c028812645a122e22ea20816678df02967c1b23bd72. DET
+    // is (𝑣₁, 𝑣₂) matrix determinant. For more details see
+    // https://www.iacr.org/archive/crypto2001/21390189.pdf
+    static constexpr auto X1 = 64502973549206556628585045361533709077_u512;
+    // Y1 should be negative, hence we calculate the determinant below adding operands instead of
+    // subtracting.
+    static constexpr auto Y1 = 303414439467246543595250775667605759171_u512;
+    static constexpr auto X2 = 367917413016453100223835821029139468248_u512;
+    static constexpr auto Y2 = 64502973549206556628585045361533709077_u512;
+    // For secp256k1 the determinant equals curve order.
+    static constexpr auto DET = uint512(Order);
+};
+// For bn254 curve and β ∈ 𝔽ₚ endomorphism ϕ : E₂ → E₂ defined as (𝑥,𝑦) → (β𝑥,𝑦) calculates [λ](𝑥,𝑦)
+// with only one multiplication in 𝔽ₚ. BETA value in Montgomery form;
+static const auto BETA =
+    Fp.to_mont(0x7ae96a2b657c07106e64479eac3434e99cf0497512f58995c1396c28719501ee_u256);
+
 }  // namespace
 
 // FIXME: Change to "uncompress_point".
@@ -117,9 +140,22 @@ std::optional<Point> secp256k1_ecdsa_recover(
     // 6. Calculate public key point Q.
     const auto R = ecc::to_proj(Fp, {r, y});
     const auto pG = ecc::to_proj(Fp, G);
-    const auto T1 = ecc::mul(Fp, pG, u1, B3);
-    const auto T2 = ecc::mul(Fp, R, u2, B3);
-    const auto pQ = ecc::add(Fp, T1, T2, B3);
+
+    const auto [u1k1, u1k2] = ecc::decompose<Config>(u1);
+    const auto [u2k1, u2k2] = ecc::decompose<Config>(u2);
+
+    const ecc::ProjPoint<uint256> pLG = {
+        Fp.mul(BETA, pG.x), !u1k2.first ? pG.y : Fp.sub(0, pG.y), pG.z};
+    const ecc::ProjPoint<uint256> pLR = {
+        Fp.mul(BETA, R.x), !u2k2.first ? R.y : Fp.sub(0, R.y), R.z};
+
+    const auto pQ = ecc::add(Fp,
+        shamir_multiply(Fp, B3, u1k1.second,
+            !u1k1.first ? pG : ecc::ProjPoint<uint256>{pG.x, Fp.sub(0, pG.y), pG.z}, u1k2.second,
+            pLG),
+        shamir_multiply(Fp, B3, u2k1.second,
+            !u2k1.first ? R : ecc::ProjPoint<uint256>{R.x, Fp.sub(0, R.y), R.z}, u2k2.second, pLR),
+        B3);
 
     const auto Q = ecc::to_affine(Fp, field_inv, pQ);