Add specialized reduction for secp256k1

This takes advantage of the modulus being near a power of 2.

Performance improvement varies by CPU, compiler, and specific
algorithm in use, but generally ranges between 10% to 20%.

src/fuzzer/mp_redc_crandall.cpp

@@ -0,0 +1,44 @@
+/*
+* (C) 2024 Jack Lloyd
+*
+* Botan is released under the Simplified BSD License (see license.txt)
+*/
+
+#include "mp_fuzzers.h"
+
+#include <botan/bigint.h>
+#include <botan/internal/loadstor.h>
+
+void fuzz(const uint8_t in[], size_t in_len) {
+   if(in_len != 8 * sizeof(word)) {
+      return;
+   }
+
+#if BOTAN_MP_WORD_BITS == 64
+   // secp256k1 modulus
+   const word C = 0x1000003d1;
+#else
+   // 128 bit prime with largest possible C
+   const word C = 0xffffffe1;
+#endif
+
+   static const Botan::BigInt refp = Botan::BigInt::power_of_2(4 * BOTAN_MP_WORD_BITS) - C;
+   static const Botan::BigInt refp2 = refp * refp;
+
+   const auto refz = Botan::BigInt::from_bytes(std::span{in, in_len});
+
+   if(refz >= refp2) {
+      return;
+   }
+
+   const auto refc = refz % refp;
+
+   std::array<word, 8> z = {};
+   for(size_t i = 0; i != 8; ++i) {
+      z[7 - i] = Botan::load_be<word>(in, i);
+   }
+
+   const auto rc = Botan::redc_crandall<word, 4, C>(z);
+
+   compare_word_vec(rc.data(), 4, refc._data(), refc.sig_words(), "Crandall reduction");
+}

src/lib/math/mp/mp_core.h

@@ -384,7 +384,7 @@ inline constexpr void bigint_monty_maybe_sub(size_t N, W z[], W x0, const W x[],
       z[i] = word_sub(x[i], p[i], &borrow);
    }
 
-   word_sub(x0, static_cast<W>(0), &borrow);
+   borrow = (x0 - borrow) > x0;
 
    CT::conditional_assign_mem(borrow, z, x, N);
 }
@@ -418,7 +418,7 @@ inline constexpr void bigint_monty_maybe_sub(W z[N], W x0, const W x[N], const W
       }
    }
 
-   word_sub(x0, static_cast<W>(0), &borrow);
+   borrow = (x0 - borrow) > x0;
 
    CT::conditional_assign_mem(borrow, z, x, N);
 }
@@ -1094,6 +1094,57 @@ void bigint_mul(word z[],
 
 void bigint_sqr(word z[], size_t z_size, const word x[], size_t x_size, size_t x_sw, word workspace[], size_t ws_size);
 
+/**
+* Return 2**B - C
+*/
+template <WordType W, size_t N, W C>
+consteval std::array<W, N> crandall_p() {
+   static_assert(C % 2 == 1);
+   std::array<W, N> P;
+   for(size_t i = 0; i != N; ++i) {
+      P[i] = WordInfo<W>::max;
+   }
+   P[0] = WordInfo<W>::max - (C - 1);
+   return P;
+}
+
+/**
+* Reduce z modulo p = 2**B - C where C is small
+*
+* z is assumed to be at most (p-1)**2
+*
+* For details on the algorithm see
+* - Handbook of Applied Cryptography, Algorithm 14.47
+* - Guide to Elliptic Curve Cryptography, Algorithm 2.54 and Note 2.55
+*
+*/
+template <WordType W, size_t N, W C>
+constexpr std::array<W, N> redc_crandall(std::span<const W, 2 * N> z) {
+   static_assert(N >= 2);
+
+   std::array<W, N> hi = {};
+
+   // hi = hi * c + lo
+
+   W carry = 0;
+   for(size_t i = 0; i != N; ++i) {
+      hi[i] = word_madd3(z[i + N], C, z[i], &carry);
+   }
+
+   // hi += carry * C
+   word carry_c[2] = {0};
+   carry_c[0] = word_madd2(carry, C, &carry_c[1]);
+
+   carry = bigint_add2_nc(hi.data(), N, carry_c, 2);
+
+   constexpr auto P = crandall_p<W, N, C>();
+
+   std::array<W, N> r = {};
+   bigint_monty_maybe_sub<N, W>(r.data(), carry, hi.data(), P.data());
+
+   return r;
+}
+
 }  // namespace Botan
 
 #endif

src/lib/math/pcurves/pcurves_secp256k1/pcurves_secp256k1.cpp

@@ -12,9 +12,33 @@ namespace Botan::PCurve {
 
 namespace {
 
-// clang-format off
 namespace secp256k1 {
 
+template <typename Params>
+class Secp256k1Rep final {
+   public:
+      static constexpr auto P = Params::P;
+      static constexpr size_t N = Params::N;
+      typedef typename Params::W W;
+
+      static_assert(WordInfo<W>::bits >= 33);
+
+      static constexpr W C = 0x1000003d1;
+
+      constexpr static std::array<W, N> one() { return std::array<W, N>{1}; }
+
+      constexpr static std::array<W, N> redc(const std::array<W, 2 * N>& z) {
+         return redc_crandall<W, N, C>(std::span{z});
+      }
+
+      constexpr static std::array<W, N> to_rep(const std::array<W, N>& x) { return x; }
+
+      constexpr static std::array<W, N> wide_to_rep(const std::array<W, 2 * N>& x) { return redc(x); }
+
+      constexpr static std::array<W, N> from_rep(const std::array<W, N>& z) { return z; }
+};
+
+// clang-format off
 class Params final : public EllipticCurveParameters<
    "FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFEFFFFFC2F",
    "0",
@@ -24,11 +48,15 @@ class Params final : public EllipticCurveParameters<
    "483ADA7726A3C4655DA4FBFC0E1108A8FD17B448A68554199C47D08FFB10D4B8"> {
 };
 
-class Curve final : public EllipticCurve<Params> {};
+// clang-format on
 
-}
+#if BOTAN_MP_WORD_BITS == 64
+class Curve final : public EllipticCurve<Params, Secp256k1Rep> {};
+#else
+class Curve final : public EllipticCurve<Params> {};
+#endif
 
-// clang-format on
+}  // namespace secp256k1
 
 }  // namespace