Merge bitcoin/bitcoin#21590: Safegcd-based modular inverses in MuHash3072

f5883286e3 Add a fuzz test for Num3072 multiplication and inversion (Pieter Wuille)
a26ce62894 Safegcd based modular inverse for Num3072 (Pieter Wuille)
91ce8cef2d Add benchmark for MuHash finalization (Pieter Wuille)

Pull request description:

  This implements a safegcd-based modular inverse for MuHash3072. It is a fairly straightforward translation of [the libsecp256k1 implementation](https://github.com/bitcoin-core/secp256k1/pull/831), with the following changes:
  * Generic for 32-bit and 64-bit
  * Specialized for the specific MuHash3072 modulus (2^3072 - 1103717).
  * A bit more C++ish
  * Far fewer sanity checks

  A benchmark is also included for MuHash3072::Finalize. The new implementation is around 100x faster on x86_64 for me (from 5.8 ms to 57 μs); for 32-bit code the factor is likely even larger.

  For more information:
    * [Original paper](https://gcd.cr.yp.to/papers.html) by Daniel J. Bernstein and Bo-Yin Yang
    * [Implementation](https://github.com/bitcoin-core/secp256k1/pull/767) for libsecp256k1 by Peter Dettman; and the [final](https://github.com/bitcoin-core/secp256k1/pull/831) version
    * [Explanation](https://github.com/bitcoin-core/secp256k1/blob/master/doc/safegcd_implementation.md) of the algorithm using Python snippets
    * [Analysis](https://github.com/sipa/safegcd-bounds) of the maximum number of iterations the algorithm needs
     * [Formal proof in Coq](https://medium.com/blockstream/a-formal-proof-of-safegcd-bounds-695e1735a348) by Russell O'Connor (for the 256-bit version of the algorithm; here we use a 3072-bit one).

ACKs for top commit:
  achow101:
    ACK f5883286e3
  TheCharlatan:
    Re-ACK f5883286e3
  dergoegge:
    tACK f5883286e3

Tree-SHA512: 275872c61d30817a82901dee93fc7153afca55c32b72a95b8768f3fd464da1b09b36f952f30e70225e766b580751cfb9b874b2feaeb73ffaa6943c8062aee19a
This commit is contained in:
Ava Chow
2025-01-27 16:50:16 -05:00
6 changed files with 514 additions and 93 deletions

View File

@@ -2,13 +2,169 @@
// Distributed under the MIT software license, see the accompanying
// file COPYING or http://www.opensource.org/licenses/mit-license.php.
#include <arith_uint256.h>
#include <crypto/muhash.h>
#include <span.h>
#include <uint256.h>
#include <test/fuzz/FuzzedDataProvider.h>
#include <test/fuzz/fuzz.h>
#include <test/fuzz/util.h>
#include <algorithm>
#include <array>
#include <vector>
namespace {
/** Class to represent 6144-bit numbers using arith_uint256 code.
*
* 6144 is sufficient to represent the product of two 3072-bit numbers. */
class arith_uint6144 : public base_uint<6144> {
public:
arith_uint6144(uint64_t x) : base_uint{x} {}
/** Construct an arith_uint6144 from any multiple of 4 bytes in LE notation,
* up to 768 bytes. */
arith_uint6144(Span<const uint8_t> bytes) : base_uint{}
{
assert(bytes.size() % 4 == 0);
assert(bytes.size() <= 768);
for (unsigned i = 0; i * 4 < bytes.size(); ++i) {
pn[i] = ReadLE32(bytes.data() + 4 * i);
}
}
/** Serialize an arithm_uint6144 to any multiply of 4 bytes in LE notation,
* on the condition that the represented number fits. */
void Serialize(Span<uint8_t> bytes) {
assert(bytes.size() % 4 == 0);
assert(bytes.size() <= 768);
for (unsigned i = 0; i * 4 < bytes.size(); ++i) {
WriteLE32(bytes.data() + 4 * i, pn[i]);
}
for (unsigned i = bytes.size() / 4; i * 4 < 768; ++i) {
assert(pn[i] == 0);
}
};
};
/** The MuHash3072 modulus (2**3072 - 1103717) as 768 LE8 bytes. */
constexpr std::array<const uint8_t, 768> MODULUS_BYTES = {
155, 40, 239, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255, 255,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0
};
const arith_uint6144 ZERO{0};
const arith_uint6144 ONE{1};
const arith_uint6144 MODULUS{MODULUS_BYTES};
/** Update value to be the modulus of the input modulo MODULUS. */
void Reduce(arith_uint6144& value)
{
arith_uint6144 tmp = value;
tmp /= MODULUS;
tmp *= MODULUS;
value -= tmp;
}
} // namespace
FUZZ_TARGET(num3072_mul)
{
// Test multiplication
FuzzedDataProvider provider{buffer.data(), buffer.size()};
// Read two 3072-bit numbers from fuzz input, and construct arith_uint6144
// and Num3072 objects with the read values.
uint16_t data_a_len = provider.ConsumeIntegralInRange(0, 384);
uint8_t data_a[384] = {0};
provider.ConsumeData(data_a, data_a_len);
arith_uint6144 a_uint{data_a};
Num3072 a_num{data_a};
uint8_t data_b[384] = {0};
provider.ConsumeData(data_b, 384);
arith_uint6144 b_uint{data_b};
Num3072 b_num{data_b};
// Multiply the first number with the second, in both representations.
a_num.Multiply(b_num);
a_uint *= b_uint;
Reduce(a_uint);
// Serialize both to bytes and compare.
uint8_t buf_num[384], buf_uint[384];
a_num.ToBytes(buf_num);
a_uint.Serialize(buf_uint);
assert(std::ranges::equal(buf_num, buf_uint));
}
FUZZ_TARGET(num3072_inv)
{
// Test inversion
FuzzedDataProvider provider{buffer.data(), buffer.size()};
// Read a 3072-bit number from fuzz input, and construct arith_uint6144
// and Num3072 objects with the read values.
uint8_t data[384] = {0};
provider.ConsumeData(data, 384);
Num3072 num{data};
arith_uint6144 uint{data};
// Bail out if the number has no inverse.
if ((uint == ZERO) || (uint == MODULUS)) return;
// Compute the inverse of the Num3072 object.
Num3072 inv;
inv.SetToOne();
inv.Divide(num);
// Convert the computed inverse to arith_uint6144.
uint8_t buf[384];
inv.ToBytes(buf);
arith_uint6144 uint_inv{buf};
// Multiply the original and the inverse, and expect 1.
uint *= uint_inv;
Reduce(uint);
assert(uint == ONE);
}
FUZZ_TARGET(muhash)
{
FuzzedDataProvider fuzzed_data_provider{buffer.data(), buffer.size()};