crypto: Pseudo-Mersenne field multiplication for secp256k1

include/evmmax/evmmax.hpp

@@ -59,6 +59,48 @@ constexpr std::pair<uint64_t, uint64_t> addmul(
     return {p[1], p[0]};
 }
 
+/// Returns the c of a "pseudo-Mersenne" modulus 2ⁿ-c, where n is the bit width of the type.
+/// Returns 0 if c does not fit in a single word, i.e. if the modulus is not pseudo-Mersenne
+/// for the purpose of pseudo_mersenne_reduce().
+template <typename UintT>
+consteval uint64_t pseudo_mersenne_c(const UintT& mod) noexcept
+{
+    const UintT c = -mod;  // 2ⁿ - mod, because the negation is modulo 2ⁿ.
+    return (c == UintT{c[0]}) ? c[0] : 0;
+}
+
+/// Reduces a double-width value modulo the pseudo-Mersenne modulus Mod = 2ⁿ-c.
+/// Accepts any value of the double width and returns the fully reduced result.
+template <const auto& Mod>
+constexpr auto pseudo_mersenne_reduce(
+    const intx::uint<2 * std::remove_cvref_t<decltype(Mod)>::num_bits>& p) noexcept
+{
+    using UintT = std::remove_cvref_t<decltype(Mod)>;
+    constexpr auto S = UintT::num_words;
+    constexpr auto C = pseudo_mersenne_c(Mod);
+    static_assert(C != 0, "the modulus is not pseudo-Mersenne");
+
+    // Fold the high half into the low one: h⋅2ⁿ + l ≡ l + h⋅c (mod 2ⁿ-c).
+    UintT t;
+    uint64_t c = 0;
+#pragma GCC unroll 8
+    for (size_t i = 0; i != S; ++i)
+        std::tie(c, t[i]) = addmul(p[i], p[S + i], C, c);
+
+    // Fold the leftover word the same way. It is not greater than c, so the product spans
+    // two words at most and the addition below overflows by at most 1.
+    auto [r, overflow] = addc(t, UintT{intx::umul(c, C)});
+
+    // Fold the overflow again. This cannot overflow, because the addition above only overflows
+    // when the low part of its result is as small as the two-word product.
+    if (overflow) [[unlikely]]
+        r += C;
+
+    if (r >= Mod) [[unlikely]]  // The result may exceed the modulus, but by less than c.
+        r -= Mod;
+    return r;
+}
+
 /// The modular arithmetic operations for EVMMAX (EVM Modular Arithmetic Extensions).
 template <typename UintT>
 class ModArith
@@ -158,9 +200,9 @@ class ModArith
         return (d.carry) ? s : d.value;
     }
 
-    /// Computes modular inverse of x in Montgomery form. Result is in Montgomery form.
+    /// Computes u⋅x⁻¹ % mod for the given initial value of the Bézout coefficient u.
     /// Returns 0 for non-invertible x (including x == 0).
-    constexpr UintT inv(const UintT& x) const noexcept
+    constexpr UintT inv_scaled(const UintT& x, UintT u) const noexcept
     {
         assert((mod_ & 1) == 1);
         assert(mod_ >= 3);
@@ -179,12 +221,6 @@ class ModArith
         // TODO: The same paper has additional optimizations that could be applied.
         UintT a = x;
         UintT b = mod_;
-
-        // Bézout's coefficients are originally initialized to 1 and 0. But because the input x
-        // is in Montgomery form XR the algorithm would compute X⁻¹R⁻¹. To get the expected X⁻¹R,
-        // we need to multiply the result by R². We can achieve the same effect "for free"
-        // by initializing u to R² instead of 1.
-        UintT u = r_squared_;
         UintT v = 0;
 
         while (a != 0)
@@ -223,5 +259,106 @@ class ModArith
             v = 0;  // not invertible
         return v;
     }
+
+    /// Computes modular inverse of x in Montgomery form. Result is in Montgomery form.
+    /// Returns 0 for non-invertible x (including x == 0).
+    constexpr UintT inv(const UintT& x) const noexcept
+    {
+        // Bézout's coefficient u is originally initialized to 1. But because the input x is in
+        // Montgomery form XR the algorithm would compute X⁻¹R⁻¹. To get the expected X⁻¹R, we
+        // need to multiply the result by R². We can achieve the same effect "for free"
+        // by initializing u to R² instead of 1.
+        return inv_scaled(x, r_squared_);
+    }
+};
+
+/// The modular arithmetic backend of a prime field.
+///
+/// An implementation keeps values in an internal representation of its choice. Every operation
+/// takes and returns fully reduced values, i.e. less than the modulus. Additionally,
+/// to_internal() requires its argument to be already reduced, because a plain representation
+/// cannot reduce it. inv() returns 0 for non-invertible input, including 0.
+template <typename T>
+concept FieldArith = requires(const typename T::uint_type& x, const typename T::uint_type& y) {
+    { T::to_internal(x) } -> std::same_as<typename T::uint_type>;
+    { T::from_internal(x) } -> std::same_as<typename T::uint_type>;
+    { T::mul(x, y) } -> std::same_as<typename T::uint_type>;
+    { T::add(x, y) } -> std::same_as<typename T::uint_type>;
+    { T::sub(x, y) } -> std::same_as<typename T::uint_type>;
+    { T::inv(x) } -> std::same_as<typename T::uint_type>;
 };
+
+/// The Montgomery multiplication backend, usable with any odd modulus.
+/// Values are kept in the Montgomery form.
+template <const auto& Mod>
+struct MontgomeryArith
+{
+    using uint_type = std::remove_cvref_t<decltype(Mod)>;
+
+    static constexpr uint_type to_internal(const uint_type& v) noexcept { return M.to_mont(v); }
+    static constexpr uint_type from_internal(const uint_type& v) noexcept { return M.from_mont(v); }
+
+    static constexpr uint_type mul(const uint_type& x, const uint_type& y) noexcept
+    {
+        return M.mul(x, y);
+    }
+    static constexpr uint_type add(const uint_type& x, const uint_type& y) noexcept
+    {
+        return M.add(x, y);
+    }
+    static constexpr uint_type sub(const uint_type& x, const uint_type& y) noexcept
+    {
+        return M.sub(x, y);
+    }
+    static constexpr uint_type inv(const uint_type& x) noexcept { return M.inv(x); }
+
+private:
+    static constexpr ModArith<uint_type> M{Mod};
+};
+
+/// The pseudo-Mersenne backend for a modulus 2ⁿ-c with a single-word c.
+/// Values are kept plain, so the conversions to and from the internal form are free.
+template <const auto& Mod>
+struct PseudoMersenneArith
+{
+    using uint_type = std::remove_cvref_t<decltype(Mod)>;
+
+    static constexpr uint_type to_internal(const uint_type& v) noexcept
+    {
+        // Unlike the conversion to the Montgomery form, this does not reduce the input.
+        assert(v < Mod);
+        return v;
+    }
+    static constexpr uint_type from_internal(const uint_type& v) noexcept { return v; }
+
+    static constexpr uint_type mul(const uint_type& x, const uint_type& y) noexcept
+    {
+        return pseudo_mersenne_reduce<Mod>(intx::umul(x, y));
+    }
+    static constexpr uint_type add(const uint_type& x, const uint_type& y) noexcept
+    {
+        return M.add(x, y);
+    }
+    static constexpr uint_type sub(const uint_type& x, const uint_type& y) noexcept
+    {
+        return M.sub(x, y);
+    }
+
+    /// The plain representation needs no R² factor folded in, unlike the Montgomery one.
+    static constexpr uint_type inv(const uint_type& x) noexcept
+    {
+        return M.inv_scaled(x, uint_type{1});
+    }
+
+private:
+    /// Reused for the operations not affected by the representation and for the binary GCD.
+    static constexpr ModArith<uint_type> M{Mod};
+};
+
+/// Selects the arithmetic backend from the structure of the modulus.
+/// Both alternatives are named for any modulus, but only the selected one is instantiated,
+/// so a backend must not place its precondition in the class scope.
+template <const auto& Mod>
+using ArithFor =
+    std::conditional_t<pseudo_mersenne_c(Mod) != 0, PseudoMersenneArith<Mod>, MontgomeryArith<Mod>>;
 }  // namespace evmmax

lib/evmone_precompiles/ecc.hpp

@@ -30,11 +30,15 @@ template <FieldSpec Spec>
 class FieldElement
 {
     using uint_type = std::remove_const_t<decltype(Spec::ORDER)>;
-    static constexpr ModArith<uint_type> Fp{Spec::ORDER};
+
+    /// The arithmetic backend, selected by the modulus. It also decides the internal
+    /// representation of the value, so this class must not assume any.
+    using Arith = ArithFor<Spec::ORDER>;
+    static_assert(FieldArith<Arith>);
 
     uint_type value_;
 
-    /// Wraps a value into the Element type assuming it is already in the internal ModArith form.
+    /// Wraps a value into the Element type assuming it is already in the internal form.
     [[gnu::always_inline]] static constexpr FieldElement wrap(const uint_type& v) noexcept
     {
         FieldElement element;
@@ -48,9 +52,9 @@ class FieldElement
 
     FieldElement() = default;
 
-    constexpr explicit FieldElement(uint_type v) : value_{Fp.to_mont(v)} {}
+    constexpr explicit FieldElement(uint_type v) : value_{Arith::to_internal(v)} {}
 
-    constexpr uint_type value() const noexcept { return Fp.from_mont(value_); }
+    constexpr uint_type value() const noexcept { return Arith::from_internal(value_); }
 
     /// The valid range for from_bytes().
     enum class Range : bool
@@ -85,45 +89,45 @@ class FieldElement
 
     friend constexpr auto operator*(const FieldElement& a, const FieldElement& b) noexcept
     {
-        return wrap(Fp.mul(a.value_, b.value_));
+        return wrap(Arith::mul(a.value_, b.value_));
     }
 
     friend constexpr auto operator+(const FieldElement& a, const FieldElement& b) noexcept
     {
-        return wrap(Fp.add(a.value_, b.value_));
+        return wrap(Arith::add(a.value_, b.value_));
     }
 
     FieldElement& operator+=(const FieldElement& b) noexcept
     {
-        value_ = Fp.add(value_, b.value_);
+        value_ = Arith::add(value_, b.value_);
         return *this;
     }
 
     friend constexpr auto operator-(const FieldElement& a, const FieldElement& b) noexcept
     {
-        return wrap(Fp.sub(a.value_, b.value_));
+        return wrap(Arith::sub(a.value_, b.value_));
     }
 
     friend constexpr auto operator-(const FieldElement& a) noexcept
     {
-        return wrap(Fp.sub(0, a.value_));
+        return wrap(Arith::sub(0, a.value_));
     }
 
-    /// Division returns 0 when the divisor is 0. See ModArith::inv().
+    /// Division returns 0 when the divisor is 0. See the FieldArith concept.
     friend constexpr auto operator/(one_t, const FieldElement& a) noexcept
     {
-        return wrap(Fp.inv(a.value_));
+        return wrap(Arith::inv(a.value_));
     }
 
-    /// Division returns 0 when the divisor is 0. See ModArith::inv().
+    /// Division returns 0 when the divisor is 0. See the FieldArith concept.
     friend constexpr auto operator/(const FieldElement& a, const FieldElement& b) noexcept
     {
-        return wrap(Fp.mul(a.value_, Fp.inv(b.value_)));
+        return wrap(Arith::mul(a.value_, Arith::inv(b.value_)));
     }
 
     /// Named 1/x inversion method. Needed in the pairing templates.
-    /// Returns 0 when this element is 0. See ModArith::inv().
-    constexpr auto inv() const noexcept { return wrap(Fp.inv(value_)); }
+    /// Returns 0 when this element is 0. See the FieldArith concept.
+    constexpr auto inv() const noexcept { return wrap(Arith::inv(value_)); }
 
     /// Named one element. Needed in the pairing templates.
     static constexpr auto one() noexcept { return FieldElement{1}; }

test/unittests/evmmax_test.cpp

@@ -5,6 +5,7 @@
 #include <evmmax/evmmax.hpp>
 #include <gtest/gtest.h>
 #include <array>
+#include <vector>
 
 using namespace intx;
 using namespace evmmax;
@@ -18,6 +19,13 @@ constexpr auto BLS12384Mod =
     0x1a0111ea397fe69a4b1ba7b6434bacd764774b84f38512bf6730d2a0f6b0f6241eabfffeb153ffffb9feffffffffaaab_u384;
 
 
+// Only the secp256k1 field prime (2²⁵⁶-2³²-977) and 2²⁵⁶-1 are pseudo-Mersenne.
+static_assert(pseudo_mersenne_c(Secp256k1Mod) == 0x1000003d1);
+static_assert(pseudo_mersenne_c(M256) == 1);
+static_assert(pseudo_mersenne_c(P23) == 0);
+static_assert(pseudo_mersenne_c(BN254Mod) == 0);
+static_assert(pseudo_mersenne_c(BLS12384Mod) == 0);
+
 template <typename UintT, const UintT& Mod>
 struct ModA : ModArith<UintT>
 {
@@ -167,3 +175,47 @@ TYPED_TEST(evmmax_test, inv)
         EXPECT_EQ(m.from_mont(pm), 1);
     }
 }
+
+/// Reduces by generic division, as the reference for pseudo_mersenne_reduce().
+template <const auto& Mod>
+static auto reduce_ref(
+    const intx::uint<2 * std::remove_cvref_t<decltype(Mod)>::num_bits>& p) noexcept
+{
+    using UintT = std::remove_cvref_t<decltype(Mod)>;
+    return static_cast<UintT>(udivrem(p, decltype(p){Mod}).rem);
+}
+
+/// Checks pseudo_mersenne_reduce() against the reference for the whole input range, including
+/// the two carry foldings that no product of two canonical operands can reach.
+template <const auto& Mod>
+static void test_reduce()
+{
+    using UintT = std::remove_cvref_t<decltype(Mod)>;
+    using WideT = intx::uint<2 * UintT::num_bits>;
+    static constexpr auto C = pseudo_mersenne_c(Mod);
+
+    std::vector<WideT> inputs{0, 1, WideT{Mod} - 1, WideT{Mod} + 1, ~WideT{},
+        umul(Mod - 1, Mod - 1), WideT{~UintT{}}, WideT{~UintT{}} - C, WideT{C}};
+
+    // Reaches the final subtraction: a zero high half makes the folding an identity, so any
+    // value in [Mod, 2ⁿ) needs the subtraction. Mod itself must reduce to 0.
+    inputs.push_back(WideT{Mod});
+
+    // Reaches the folding of the leftover word's overflow: solving l + h⋅C = 2ⁿ + (2ⁿ-1) leaves
+    // the first folding at 2ⁿ-1 with a leftover word of 1, so adding 1⋅C overflows.
+    const auto target = (WideT{1} << (UintT::num_bits + 1)) - 1;
+    inputs.push_back(((target / C) << UintT::num_bits) | (target % C));
+
+    for (const auto& p : inputs)
+        EXPECT_EQ(pseudo_mersenne_reduce<Mod>(p), reduce_ref<Mod>(p)) << to_string(p, 16);
+}
+
+TEST(evmmax, pseudo_mersenne_reduce_secp256k1)
+{
+    test_reduce<Secp256k1Mod>();
+}
+
+TEST(evmmax, pseudo_mersenne_reduce_m256)
+{
+    test_reduce<M256>();
+}