evmmax: Implement ate pairing for bn254 curve (#852)

Implement optimized version of ecparing precompile in evmone according
to [EIP197](https://eips.ethereum.org/EIPS/eip-197) and based on
[High-Speed Software Implementation of the Optimal Ate Pairing over
Barreto–Naehrig Curves](https://eprint.iacr.org/2010/354.pdf)

- Extend the interface with Jacobian caordinates point 
- Add template struct which implements basic logic for base and
extension fields.
- Add helper functions needed by the ate pairing.
- Implement miller loop and final exponentiation.
- Add unit tests and benchmark.
- Enable the implementation to be used by state and blockchain tests

Author: chfast

lib/evmone_precompiles/CMakeLists.txt

@@ -16,6 +16,10 @@ target_sources(
     bn254.hpp
     bn254.cpp
     ecc.hpp
+    pairing/bn254/fields.hpp
+    pairing/bn254/pairing.cpp
+    pairing/bn254/utils.hpp
+    pairing/field_template.hpp
     ripemd160.hpp
     ripemd160.cpp
     secp256k1.hpp

lib/evmone_precompiles/bn254.hpp

@@ -4,6 +4,9 @@
 #pragma once
 
 #include "ecc.hpp"
+#include <optional>
+#include <span>
+#include <vector>
 
 namespace evmmax::bn254
 {
@@ -14,6 +17,7 @@ inline constexpr auto FieldPrime =
     0x30644e72e131a029b85045b68181585d97816a916871ca8d3c208c16d87cfd47_u256;
 
 using Point = ecc::Point<uint256>;
+using ExtPoint = ecc::Point<std::pair<uint256, uint256>>;
 
 /// Validates that point is from the bn254 curve group
 ///
@@ -36,4 +40,12 @@ Point add(const Point& pt1, const Point& pt2) noexcept;
 /// Computes [c]P for a point in affine coordinate on the bn254 curve,
 Point mul(const Point& pt, const uint256& c) noexcept;
 
+/// ate paring implementation for bn254 curve according to https://eips.ethereum.org/EIPS/eip-197
+///
+/// @param pairs  Sequence of point pairs: a point from the bn254 curve G1 group over the base field
+///               followed by a point from twisted curve G2 group over extension field Fq^2.
+/// @return       `true` when  ∏e(vG2[i], vG1[i]) == 1 for i in [0, n] else `false`.
+///               std::nullopt on error.
+std::optional<bool> pairing_check(std::span<const std::pair<Point, ExtPoint>> pairs) noexcept;
+
 }  // namespace evmmax::bn254

lib/evmone_precompiles/ecc.hpp

@@ -9,14 +9,16 @@ namespace evmmax::ecc
 {
 
 /// The affine (two coordinates) point on an Elliptic Curve over a prime field.
-template <typename IntT>
+template <typename ValueT>
 struct Point
 {
-    IntT x = 0;
-    IntT y = 0;
+    ValueT x = {};
+    ValueT y = {};
 
     friend constexpr bool operator==(const Point& a, const Point& b) noexcept = default;
 
+    friend constexpr Point operator-(const Point& p) noexcept { return {p.x, -p.y}; }
+
     /// Checks if the point represents the special "infinity" value.
     [[nodiscard]] constexpr bool is_inf() const noexcept { return *this == Point{}; }
 };
@@ -32,10 +34,41 @@ struct ProjPoint
 
     /// Checks if the point represents the special "infinity" value.
     [[nodiscard]] constexpr bool is_inf() const noexcept { return x == 0 && z == 0; }
+
+    friend constexpr ProjPoint operator-(const ProjPoint& p) noexcept { return {p.x, -p.y, p.z}; }
 };
 
 static_assert(ProjPoint<unsigned>{}.is_inf());
 
+// Jacobian (three) coordinates point implementation.
+template <typename ValueT>
+struct JacPoint
+{
+    ValueT x = 1;
+    ValueT y = 1;
+    ValueT z = 0;
+
+    // Compares two Jacobian coordinates points
+    friend constexpr bool operator==(const JacPoint& a, const JacPoint& b) noexcept
+    {
+        const auto bz2 = b.z * b.z;
+        const auto az2 = a.z * a.z;
+
+        const auto bz3 = bz2 * b.z;
+        const auto az3 = az2 * a.z;
+
+        return a.x * bz2 == b.x * az2 && a.y * bz3 == b.y * az3;
+    }
+
+    friend constexpr JacPoint operator-(const JacPoint& p) noexcept { return {p.x, -p.y, p.z}; }
+
+    // Creates Jacobian coordinates point from affine point
+    static constexpr JacPoint from(const ecc::Point<ValueT>& ap) noexcept
+    {
+        return {ap.x, ap.y, ValueT::one()};
+    }
+};
+
 template <typename IntT>
 using InvFn = IntT (*)(const ModArith<IntT>&, const IntT& x) noexcept;
 

lib/evmone_precompiles/pairing/bn254/fields.hpp

@@ -0,0 +1,182 @@
+// evmone: Fast Ethereum Virtual Machine implementation
+// Copyright 2024 The evmone Authors.
+// SPDX-License-Identifier: Apache-2.0
+#pragma once
+
+#include "../../bn254.hpp"
+#include "../../ecc.hpp"
+#include "../field_template.hpp"
+
+namespace evmmax::bn254
+{
+using namespace intx;
+
+/// Specifies base field value type and modular arithmetic for bn254 curve.
+struct BaseFieldConfig
+{
+    using ValueT = uint256;
+    static constexpr auto MOD_ARITH = ModArith{FieldPrime};
+    static constexpr uint256 ONE = MOD_ARITH.to_mont(1);
+};
+using Fq = ecc::BaseFieldElem<BaseFieldConfig>;
+
+// Extension fields implemented based on https://hackmd.io/@jpw/bn254#Field-extension-towers
+
+/// Specifies Fq^2 extension field for bn254 curve. Base field extended with irreducible `u^2 + 1`
+/// polynomial over the base field. `u` is the Fq^2 element.
+struct Fq2Config
+{
+    using BaseFieldT = Fq;
+    using ValueT = Fq;
+    static constexpr auto DEGREE = 2;
+};
+using Fq2 = ecc::ExtFieldElem<Fq2Config>;
+
+/// Specifies Fq^6 extension field for bn254 curve. Fq^2 field extended with irreducible
+/// `v^3 - (9 + u)` polynomial over the Fq^2 field. `v` is the Fq^6 field element.
+struct Fq6Config
+{
+    using BaseFieldT = Fq;
+    using ValueT = Fq2;
+    static constexpr uint8_t DEGREE = 3;
+    static constexpr auto ksi = Fq2({Fq::from_int(9_u256), Fq::from_int(1_u256)});
+    static constexpr auto _3_ksi_inv = Fq2({
+        Fq::from_int(0x2b149d40ceb8aaae81be18991be06ac3b5b4c5e559dbefa33267e6dc24a138e5_u256),
+        Fq::from_int(0x9713b03af0fed4cd2cafadeed8fdf4a74fa084e52d1852e4a2bd0685c315d2_u256),
+    });
+};
+using Fq6 = ecc::ExtFieldElem<Fq6Config>;
+
+/// Specifies Fq^12 extension field for bn254 curve. Fq^6 field extended with irreducible
+/// `w^2 - v` polynomial over the Fq^2 field. `v` is the Fq^6 field element.
+/// `w` is Fq^12 field element.
+struct Fq12Config
+{
+    using BaseFieldT = Fq;
+    using ValueT = Fq6;
+
+    static constexpr uint8_t DEGREE = 2;
+};
+using Fq12 = ecc::ExtFieldElem<Fq12Config>;
+
+/// Multiplies two Fq^2 field elements
+constexpr Fq2 multiply(const Fq2& a, const Fq2& b)
+{
+    return Fq2({
+        a.coeffs[0] * b.coeffs[0] - a.coeffs[1] * b.coeffs[1],
+        a.coeffs[1] * b.coeffs[0] + a.coeffs[0] * b.coeffs[1],
+    });
+}
+
+/// Multiplies two Fq^6 field elements
+constexpr Fq6 multiply(const Fq6& a, const Fq6& b)
+{
+    const auto& a0 = a.coeffs[0];
+    const auto& a1 = a.coeffs[1];
+    const auto& a2 = a.coeffs[2];
+    const auto& b0 = b.coeffs[0];
+    const auto& b1 = b.coeffs[1];
+    const auto& b2 = b.coeffs[2];
+
+    const Fq2& ksi = Fq6Config::ksi;
+
+    const auto t0 = a0 * b0;
+    const auto t1 = a1 * b1;
+    const auto t2 = a2 * b2;
+
+    const auto c0 = ((a1 + a2) * (b1 + b2) - t1 - t2) * ksi + t0;
+    const auto c1 = (a0 + a1) * (b0 + b1) - t0 - t1 + ksi * t2;
+    const auto c2 = (a0 + a2) * (b0 + b2) - t0 - t2 + t1;
+
+    return Fq6({c0, c1, c2});
+}
+
+/// Multiplies two Fq^12 field elements
+constexpr Fq12 multiply(const Fq12& a, const Fq12& b)
+{
+    const auto& a0 = a.coeffs[0];
+    const auto& a1 = a.coeffs[1];
+    const auto& b0 = b.coeffs[0];
+    const auto& b1 = b.coeffs[1];
+
+    const auto t0 = a0 * b0;
+    const auto t1 = a1 * b1;
+
+    const Fq2& ksi = Fq6Config::ksi;
+
+    const auto c0 = t0 + Fq6({ksi * t1.coeffs[2], t1.coeffs[0], t1.coeffs[1]});  // gamma is sparse.
+    const auto c1 = (a0 + a1) * (b0 + b1) - t0 - t1;
+
+    return Fq12({c0, c1});
+}
+
+/// Inverses the base field element
+inline Fq inverse(const Fq& x)
+{
+    return Fq(field_inv(BaseFieldConfig::MOD_ARITH, x.value()));
+}
+
+/// Inverses the Fq^2 field element
+inline Fq2 inverse(const Fq2& f)
+{
+    const auto& a0 = f.coeffs[0];
+    const auto& a1 = f.coeffs[1];
+    auto t0 = a0 * a0;
+    auto t1 = a1 * a1;
+
+    t0 = t0 + t1;
+    t1 = t0.inv();
+
+    const auto c0 = a0 * t1;
+    const auto c1 = -(a1 * t1);
+
+    return Fq2({c0, c1});
+}
+
+/// Inverses the Fq^6 field element
+inline Fq6 inverse(const Fq6& f)
+{
+    const auto& a0 = f.coeffs[0];
+    const auto& a1 = f.coeffs[1];
+    const auto& a2 = f.coeffs[2];
+
+    const Fq2& ksi = Fq6Config::ksi;
+
+    const auto t0 = a0 * a0;
+    const auto t1 = a1 * a1;
+    const auto t2 = a2 * a2;
+
+    const auto t3 = a0 * a1;
+    const auto t4 = a0 * a2;
+    const auto t5 = a2 * a1;
+
+    const auto c0 = t0 - ksi * t5;
+    const auto c1 = ksi * t2 - t3;
+    const auto c2 = t1 - t4;
+
+    const auto t = a0 * c0 + (a2 * c1 + a1 * c2) * ksi;
+    const auto t6 = t.inv();
+
+    return Fq6({c0 * t6, c1 * t6, c2 * t6});
+}
+
+/// Inverses the Fq^12 field element
+inline Fq12 inverse(const Fq12& f)
+{
+    const auto& a0 = f.coeffs[0];
+    const auto& a1 = f.coeffs[1];
+
+    auto t0 = a0 * a0;
+    auto t1 = a1 * a1;
+
+    const Fq2& ksi = Fq6Config::ksi;
+
+    t0 = t0 - Fq6({ksi * t1.coeffs[2], t1.coeffs[0], t1.coeffs[1]});  // gamma is sparse.
+    t1 = t0.inv();
+
+    const auto c0 = a0 * t1;
+    const auto c1 = -(a1 * t1);
+
+    return Fq12({c0, c1});
+}
+}  // namespace evmmax::bn254