diff --git a/include/evmmax/evmmax.hpp b/include/evmmax/evmmax.hpp index 9658100702..d67c04ec89 100644 --- a/include/evmmax/evmmax.hpp +++ b/include/evmmax/evmmax.hpp @@ -59,6 +59,48 @@ constexpr std::pair 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 +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 +constexpr auto pseudo_mersenne_reduce( + const intx::uint<2 * std::remove_cvref_t::num_bits>& p) noexcept +{ + using UintT = std::remove_cvref_t; + 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 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 +concept FieldArith = requires(const typename T::uint_type& x, const typename T::uint_type& y) { + { T::to_internal(x) } -> std::same_as; + { T::from_internal(x) } -> std::same_as; + { T::mul(x, y) } -> std::same_as; + { T::add(x, y) } -> std::same_as; + { T::sub(x, y) } -> std::same_as; + { T::inv(x) } -> std::same_as; }; + +/// The Montgomery multiplication backend, usable with any odd modulus. +/// Values are kept in the Montgomery form. +template +struct MontgomeryArith +{ + using uint_type = std::remove_cvref_t; + + 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 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 +struct PseudoMersenneArith +{ + using uint_type = std::remove_cvref_t; + + 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(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 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 +using ArithFor = + std::conditional_t, MontgomeryArith>; } // namespace evmmax diff --git a/lib/evmone_precompiles/ecc.hpp b/lib/evmone_precompiles/ecc.hpp index bf95cc9752..b5bfdda9d5 100644 --- a/lib/evmone_precompiles/ecc.hpp +++ b/lib/evmone_precompiles/ecc.hpp @@ -30,11 +30,15 @@ template class FieldElement { using uint_type = std::remove_const_t; - static constexpr ModArith 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; + static_assert(FieldArith); 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}; } diff --git a/test/unittests/evmmax_test.cpp b/test/unittests/evmmax_test.cpp index a721e0d7d4..80a5567d5f 100644 --- a/test/unittests/evmmax_test.cpp +++ b/test/unittests/evmmax_test.cpp @@ -5,6 +5,7 @@ #include #include #include +#include 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 struct ModA : ModArith { @@ -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 +static auto reduce_ref( + const intx::uint<2 * std::remove_cvref_t::num_bits>& p) noexcept +{ + using UintT = std::remove_cvref_t; + return static_cast(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 +static void test_reduce() +{ + using UintT = std::remove_cvref_t; + using WideT = intx::uint<2 * UintT::num_bits>; + static constexpr auto C = pseudo_mersenne_c(Mod); + + std::vector 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(p), reduce_ref(p)) << to_string(p, 16); +} + +TEST(evmmax, pseudo_mersenne_reduce_secp256k1) +{ + test_reduce(); +} + +TEST(evmmax, pseudo_mersenne_reduce_m256) +{ + test_reduce(); +}