Merge #19055: Add MuHash3072 implementation

9815332d5158d69a94abeaf465a2c07bd8e43359 test: Change MuHash Python implementation to match cpp version again (Fabian Jahr)
01297fb3ca57e4b8cbc5a89fc7c6367de33b0bc6 fuzz: Add MuHash consistency fuzz test (Fabian Jahr)
b111410914041b72961536c3e4037eba103a8085 test: Add MuHash3072 fuzz test (Fabian Jahr)
c1225273857f9fa2e2276396e3f8b3ea48306df3 bench: Add Muhash benchmarks (Fabian Jahr)
7b1242229d1fcc9277238a3aefb3431061c82bfa test: Add MuHash3072 unit tests (Fabian Jahr)
adc708c98dbf03b1735edc91f813a36580781a95 crypto: Add MuHash3072 implementation (Fabian Jahr)
0b4d290bf5b0a4d156c523431bf89aaa9ffe92e5 crypto: Add Num3072 implementation (Fabian Jahr)
589f958662a2dcaacdb9a66f1088c74828a39577 build: Check for 128 bit integer support (Fabian Jahr)

Pull request description:

  This is the first split of #18000 which implements the Muhash algorithm and uses it to calculate the UTXO set hash in `gettxoutsetinfo`.

ACKs for top commit:
  laanwj:
    Code review ACK 9815332d5158d69a94abeaf465a2c07bd8e43359

Tree-SHA512: 4bc090738f0e3d80b74bdd8122e24a8ce80121120fd37c7e4335a73e7ba4fcd7643f2a2d559e2eebf54b8e3a3bd5f12cfb27ba61ded135fda210a07a233eae45

8 files changed, 683 insertions, 1 deletions

diff --git a/src/Makefile.am b/src/Makefile.am
index 23d790d552..02797b19bc 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -408,6 +408,8 @@ crypto_libbitcoin_crypto_base_a_SOURCES = \
   crypto/hmac_sha512.h \
   crypto/poly1305.h \
   crypto/poly1305.cpp \
+  crypto/muhash.h \
+  crypto/muhash.cpp \
   crypto/ripemd160.cpp \
   crypto/ripemd160.h \
   crypto/sha1.cpp \
diff --git a/src/Makefile.test.include b/src/Makefile.test.include
index 3faa5ac968..bdaddb4f0d 100644
--- a/src/Makefile.test.include
+++ b/src/Makefile.test.include
@@ -230,6 +230,7 @@ test_fuzz_fuzz_SOURCES = \
  test/fuzz/locale.cpp \
  test/fuzz/merkleblock.cpp \
  test/fuzz/message.cpp \
+ test/fuzz/muhash.cpp \
  test/fuzz/multiplication_overflow.cpp \
  test/fuzz/net.cpp \
  test/fuzz/net_permissions.cpp \
diff --git a/src/bench/crypto_hash.cpp b/src/bench/crypto_hash.cpp
index 65d16d47d8..30fe11be6b 100644
--- a/src/bench/crypto_hash.cpp
+++ b/src/bench/crypto_hash.cpp
@@ -4,6 +4,7 @@
 
 
 #include <bench/bench.h>
+#include <crypto/muhash.h>
 #include <crypto/ripemd160.h>
 #include <crypto/sha1.h>
 #include <crypto/sha256.h>
@@ -105,6 +106,54 @@ static void FastRandom_1bit(benchmark::Bench& bench)
     });
 }
 
+static void MuHash(benchmark::Bench& bench)
+{
+    MuHash3072 acc;
+    unsigned char key[32] = {0};
+    int i = 0;
+    bench.run([&] {
+        key[0] = ++i;
+        acc *= MuHash3072(key);
+    });
+}
+
+static void MuHashMul(benchmark::Bench& bench)
+{
+    MuHash3072 acc;
+    FastRandomContext rng(true);
+    MuHash3072 muhash{rng.randbytes(32)};
+
+    bench.run([&] {
+        acc *= muhash;
+    });
+}
+
+static void MuHashDiv(benchmark::Bench& bench)
+{
+    MuHash3072 acc;
+    FastRandomContext rng(true);
+    MuHash3072 muhash{rng.randbytes(32)};
+
+    for (size_t i = 0; i < bench.epochIterations(); ++i) {
+        acc *= muhash;
+    }
+
+    bench.run([&] {
+        acc /= muhash;
+    });
+}
+
+static void MuHashPrecompute(benchmark::Bench& bench)
+{
+    MuHash3072 acc;
+    FastRandomContext rng(true);
+    std::vector<unsigned char> key{rng.randbytes(32)};
+
+    bench.run([&] {
+        MuHash3072{key};
+    });
+}
+
 BENCHMARK(RIPEMD160);
 BENCHMARK(SHA1);
 BENCHMARK(SHA256);
@@ -116,3 +165,8 @@ BENCHMARK(SipHash_32b);
 BENCHMARK(SHA256D64_1024);
 BENCHMARK(FastRandom_32bit);
 BENCHMARK(FastRandom_1bit);
+
+BENCHMARK(MuHash);
+BENCHMARK(MuHashMul);
+BENCHMARK(MuHashDiv);
+BENCHMARK(MuHashPrecompute);
diff --git a/src/crypto/muhash.cpp b/src/crypto/muhash.cpp
new file mode 100644
index 0000000000..fbd14f9325
--- /dev/null
+++ b/src/crypto/muhash.cpp
@@ -0,0 +1,338 @@
+// Copyright (c) 2017-2020 The Bitcoin Core developers
+// Distributed under the MIT software license, see the accompanying
+// file COPYING or http://www.opensource.org/licenses/mit-license.php.
+
+#include <crypto/muhash.h>
+
+#include <crypto/chacha20.h>
+#include <crypto/common.h>
+#include <hash.h>
+
+#include <cassert>
+#include <cstdio>
+#include <limits>
+
+namespace {
+
+using limb_t = Num3072::limb_t;
+using double_limb_t = Num3072::double_limb_t;
+constexpr int LIMB_SIZE = Num3072::LIMB_SIZE;
+constexpr int LIMBS = Num3072::LIMBS;
+/** 2^3072 - 1103717, the largest 3072-bit safe prime number, is used as the modulus. */
+constexpr limb_t MAX_PRIME_DIFF = 1103717;
+
+/** Extract the lowest limb of [c0,c1,c2] into n, and left shift the number by 1 limb. */
+inline void extract3(limb_t& c0, limb_t& c1, limb_t& c2, limb_t& n)
+{
+    n = c0;
+    c0 = c1;
+    c1 = c2;
+    c2 = 0;
+}
+
+/** [c0,c1] = a * b */
+inline void mul(limb_t& c0, limb_t& c1, const limb_t& a, const limb_t& b)
+{
+    double_limb_t t = (double_limb_t)a * b;
+    c1 = t >> LIMB_SIZE;
+    c0 = t;
+}
+
+/* [c0,c1,c2] += n * [d0,d1,d2]. c2 is 0 initially */
+inline void mulnadd3(limb_t& c0, limb_t& c1, limb_t& c2, limb_t& d0, limb_t& d1, limb_t& d2, const limb_t& n)
+{
+    double_limb_t t = (double_limb_t)d0 * n + c0;
+    c0 = t;
+    t >>= LIMB_SIZE;
+    t += (double_limb_t)d1 * n + c1;
+    c1 = t;
+    t >>= LIMB_SIZE;
+    c2 = t + d2 * n;
+}
+
+/* [c0,c1] *= n */
+inline void muln2(limb_t& c0, limb_t& c1, const limb_t& n)
+{
+    double_limb_t t = (double_limb_t)c0 * n;
+    c0 = t;
+    t >>= LIMB_SIZE;
+    t += (double_limb_t)c1 * n;
+    c1 = t;
+}
+
+/** [c0,c1,c2] += a * b */
+inline void muladd3(limb_t& c0, limb_t& c1, limb_t& c2, const limb_t& a, const limb_t& b)
+{
+    double_limb_t t = (double_limb_t)a * b;
+    limb_t th = t >> LIMB_SIZE;
+    limb_t tl = t;
+
+    c0 += tl;
+    th += (c0 < tl) ? 1 : 0;
+    c1 += th;
+    c2 += (c1 < th) ? 1 : 0;
+}
+
+/** [c0,c1,c2] += 2 * a * b */
+inline void muldbladd3(limb_t& c0, limb_t& c1, limb_t& c2, const limb_t& a, const limb_t& b)
+{
+    double_limb_t t = (double_limb_t)a * b;
+    limb_t th = t >> LIMB_SIZE;
+    limb_t tl = t;
+
+    c0 += tl;
+    limb_t tt = th + ((c0 < tl) ? 1 : 0);
+    c1 += tt;
+    c2 += (c1 < tt) ? 1 : 0;
+    c0 += tl;
+    th += (c0 < tl) ? 1 : 0;
+    c1 += th;
+    c2 += (c1 < th) ? 1 : 0;
+}
+
+/**
+ * Add limb a to [c0,c1]: [c0,c1] += a. Then extract the lowest
+ * limb of [c0,c1] into n, and left shift the number by 1 limb.
+ * */
+inline void addnextract2(limb_t& c0, limb_t& c1, const limb_t& a, limb_t& n)
+{
+    limb_t c2 = 0;
+
+    // add
+    c0 += a;
+    if (c0 < a) {
+        c1 += 1;
+
+        // Handle case when c1 has overflown
+        if (c1 == 0)
+            c2 = 1;
+    }
+
+    // extract
+    n = c0;
+    c0 = c1;
+    c1 = c2;
+}
+
+/** in_out = in_out^(2^sq) * mul */
+inline void square_n_mul(Num3072& in_out, const int sq, const Num3072& mul)
+{
+    for (int j = 0; j < sq; ++j) in_out.Square();
+    in_out.Multiply(mul);
+}
+
+} // namespace
+
+/** Indicates wether d is larger than the modulus. */
+bool Num3072::IsOverflow() const
+{
+    if (this->limbs[0] <= std::numeric_limits<limb_t>::max() - MAX_PRIME_DIFF) return false;
+    for (int i = 1; i < LIMBS; ++i) {
+        if (this->limbs[i] != std::numeric_limits<limb_t>::max()) return false;
+    }
+    return true;
+}
+
+void Num3072::FullReduce()
+{
+    limb_t c0 = MAX_PRIME_DIFF;
+    limb_t c1 = 0;
+    for (int i = 0; i < LIMBS; ++i) {
+        addnextract2(c0, c1, this->limbs[i], this->limbs[i]);
+    }
+}
+
+Num3072 Num3072::GetInverse() const
+{
+    // For fast exponentiation a sliding window exponentiation with repunit
+    // precomputation is utilized. See "Fast Point Decompression for Standard
+    // Elliptic Curves" (Brumley, Järvinen, 2008).
+
+    Num3072 p[12]; // p[i] = a^(2^(2^i)-1)
+    Num3072 out;
+
+    p[0] = *this;
+
+    for (int i = 0; i < 11; ++i) {
+        p[i + 1] = p[i];
+        for (int j = 0; j < (1 << i); ++j) p[i + 1].Square();
+        p[i + 1].Multiply(p[i]);
+    }
+
+    out = p[11];
+
+    square_n_mul(out, 512, p[9]);
+    square_n_mul(out, 256, p[8]);
+    square_n_mul(out, 128, p[7]);
+    square_n_mul(out, 64, p[6]);
+    square_n_mul(out, 32, p[5]);
+    square_n_mul(out, 8, p[3]);
+    square_n_mul(out, 2, p[1]);
+    square_n_mul(out, 1, p[0]);
+    square_n_mul(out, 5, p[2]);
+    square_n_mul(out, 3, p[0]);
+    square_n_mul(out, 2, p[0]);
+    square_n_mul(out, 4, p[0]);
+    square_n_mul(out, 4, p[1]);
+    square_n_mul(out, 3, p[0]);
+
+    return out;
+}
+
+void Num3072::Multiply(const Num3072& a)
+{
+    limb_t c0 = 0, c1 = 0, c2 = 0;
+    Num3072 tmp;
+
+    /* Compute limbs 0..N-2 of this*a into tmp, including one reduction. */
+    for (int j = 0; j < LIMBS - 1; ++j) {
+        limb_t d0 = 0, d1 = 0, d2 = 0;
+        mul(d0, d1, this->limbs[1 + j], a.limbs[LIMBS + j - (1 + j)]);
+        for (int i = 2 + j; i < LIMBS; ++i) muladd3(d0, d1, d2, this->limbs[i], a.limbs[LIMBS + j - i]);
+        mulnadd3(c0, c1, c2, d0, d1, d2, MAX_PRIME_DIFF);
+        for (int i = 0; i < j + 1; ++i) muladd3(c0, c1, c2, this->limbs[i], a.limbs[j - i]);
+        extract3(c0, c1, c2, tmp.limbs[j]);
+    }
+
+    /* Compute limb N-1 of a*b into tmp. */
+    assert(c2 == 0);
+    for (int i = 0; i < LIMBS; ++i) muladd3(c0, c1, c2, this->limbs[i], a.limbs[LIMBS - 1 - i]);
+    extract3(c0, c1, c2, tmp.limbs[LIMBS - 1]);
+
+    /* Perform a second reduction. */
+    muln2(c0, c1, MAX_PRIME_DIFF);
+    for (int j = 0; j < LIMBS; ++j) {
+        addnextract2(c0, c1, tmp.limbs[j], this->limbs[j]);
+    }
+
+    assert(c1 == 0);
+    assert(c0 == 0 || c0 == 1);
+
+    /* Perform up to two more reductions if the internal state has already
+     * overflown the MAX of Num3072 or if it is larger than the modulus or
+     * if both are the case.
+     * */
+    if (this->IsOverflow()) this->FullReduce();
+    if (c0) this->FullReduce();
+}
+
+void Num3072::Square()
+{
+    limb_t c0 = 0, c1 = 0, c2 = 0;
+    Num3072 tmp;
+
+    /* Compute limbs 0..N-2 of this*this into tmp, including one reduction. */
+    for (int j = 0; j < LIMBS - 1; ++j) {
+        limb_t d0 = 0, d1 = 0, d2 = 0;
+        for (int i = 0; i < (LIMBS - 1 - j) / 2; ++i) muldbladd3(d0, d1, d2, this->limbs[i + j + 1], this->limbs[LIMBS - 1 - i]);
+        if ((j + 1) & 1) muladd3(d0, d1, d2, this->limbs[(LIMBS - 1 - j) / 2 + j + 1], this->limbs[LIMBS - 1 - (LIMBS - 1 - j) / 2]);
+        mulnadd3(c0, c1, c2, d0, d1, d2, MAX_PRIME_DIFF);
+        for (int i = 0; i < (j + 1) / 2; ++i) muldbladd3(c0, c1, c2, this->limbs[i], this->limbs[j - i]);
+        if ((j + 1) & 1) muladd3(c0, c1, c2, this->limbs[(j + 1) / 2], this->limbs[j - (j + 1) / 2]);
+        extract3(c0, c1, c2, tmp.limbs[j]);
+    }
+
+    assert(c2 == 0);
+    for (int i = 0; i < LIMBS / 2; ++i) muldbladd3(c0, c1, c2, this->limbs[i], this->limbs[LIMBS - 1 - i]);
+    extract3(c0, c1, c2, tmp.limbs[LIMBS - 1]);
+
+    /* Perform a second reduction. */
+    muln2(c0, c1, MAX_PRIME_DIFF);
+    for (int j = 0; j < LIMBS; ++j) {
+        addnextract2(c0, c1, tmp.limbs[j], this->limbs[j]);
+    }
+
+    assert(c1 == 0);
+    assert(c0 == 0 || c0 == 1);
+
+    /* Perform up to two more reductions if the internal state has already
+     * overflown the MAX of Num3072 or if it is larger than the modulus or
+     * if both are the case.
+     * */
+    if (this->IsOverflow()) this->FullReduce();
+    if (c0) this->FullReduce();
+}
+
+void Num3072::SetToOne()
+{
+    this->limbs[0] = 1;
+    for (int i = 1; i < LIMBS; ++i) this->limbs[i] = 0;
+}
+
+void Num3072::Divide(const Num3072& a)
+{
+    if (this->IsOverflow()) this->FullReduce();
+
+    Num3072 inv{};
+    if (a.IsOverflow()) {
+        Num3072 b = a;
+        b.FullReduce();
+        inv = b.GetInverse();
+    } else {
+        inv = a.GetInverse();
+    }
+
+    this->Multiply(inv);
+    if (this->IsOverflow()) this->FullReduce();
+}
+
+Num3072 MuHash3072::ToNum3072(Span<const unsigned char> in) {
+    Num3072 out{};
+    uint256 hashed_in = (CHashWriter(SER_DISK, 0) << in).GetSHA256();
+    unsigned char tmp[BYTE_SIZE];
+    ChaCha20(hashed_in.data(), hashed_in.size()).Keystream(tmp, BYTE_SIZE);
+    for (int i = 0; i < LIMBS; ++i) {
+        if (sizeof(limb_t) == 4) {
+            out.limbs[i] = ReadLE32(tmp + 4 * i);
+        } else if (sizeof(limb_t) == 8) {
+            out.limbs[i] = ReadLE64(tmp + 8 * i);
+        }
+    }
+    return out;
+}
+
+MuHash3072::MuHash3072(Span<const unsigned char> in) noexcept
+{
+    m_numerator = ToNum3072(in);
+}
+
+void MuHash3072::Finalize(uint256& out) noexcept
+{
+    m_numerator.Divide(m_denominator);
+    m_denominator.SetToOne();  // Needed to keep the MuHash object valid
+
+    unsigned char data[384];
+    for (int i = 0; i < LIMBS; ++i) {
+        if (sizeof(limb_t) == 4) {
+            WriteLE32(data + i * 4, m_numerator.limbs[i]);
+        } else if (sizeof(limb_t) == 8) {
+            WriteLE64(data + i * 8, m_numerator.limbs[i]);
+        }
+    }
+
+    out = (CHashWriter(SER_DISK, 0) << data).GetSHA256();
+}
+
+MuHash3072& MuHash3072::operator*=(const MuHash3072& mul) noexcept
+{
+    m_numerator.Multiply(mul.m_numerator);
+    m_denominator.Multiply(mul.m_denominator);
+    return *this;
+}
+
+MuHash3072& MuHash3072::operator/=(const MuHash3072& div) noexcept
+{
+    m_numerator.Multiply(div.m_denominator);
+    m_denominator.Multiply(div.m_numerator);
+    return *this;
+}
+
+MuHash3072& MuHash3072::Insert(Span<const unsigned char> in) noexcept {
+    m_numerator.Multiply(ToNum3072(in));
+    return *this;
+}
+
+MuHash3072& MuHash3072::Remove(Span<const unsigned char> in) noexcept {
+    m_numerator.Divide(ToNum3072(in));
+    return *this;
+}
diff --git a/src/crypto/muhash.h b/src/crypto/muhash.h
new file mode 100644
index 0000000000..0c710007c4
--- /dev/null
+++ b/src/crypto/muhash.h
@@ -0,0 +1,130 @@
+// Copyright (c) 2017-2020 The Bitcoin Core developers
+// Distributed under the MIT software license, see the accompanying
+// file COPYING or http://www.opensource.org/licenses/mit-license.php.
+
+#ifndef BITCOIN_CRYPTO_MUHASH_H
+#define BITCOIN_CRYPTO_MUHASH_H
+
+#if defined(HAVE_CONFIG_H)
+#include <config/bitcoin-config.h>
+#endif
+
+#include <serialize.h>
+#include <uint256.h>
+
+#include <stdint.h>
+
+class Num3072
+{
+private:
+    void FullReduce();
+    bool IsOverflow() const;
+    Num3072 GetInverse() const;
+
+public:
+
+#ifdef HAVE___INT128
+    typedef unsigned __int128 double_limb_t;
+    typedef uint64_t limb_t;
+    static constexpr int LIMBS = 48;
+    static constexpr int LIMB_SIZE = 64;
+#else
+    typedef uint64_t double_limb_t;
+    typedef uint32_t limb_t;
+    static constexpr int LIMBS = 96;
+    static constexpr int LIMB_SIZE = 32;
+#endif
+    limb_t limbs[LIMBS];
+
+    // Sanity check for Num3072 constants
+    static_assert(LIMB_SIZE * LIMBS == 3072, "Num3072 isn't 3072 bits");
+    static_assert(sizeof(double_limb_t) == sizeof(limb_t) * 2, "bad size for double_limb_t");
+    static_assert(sizeof(limb_t) * 8 == LIMB_SIZE, "LIMB_SIZE is incorrect");
+
+    // Hard coded values in MuHash3072 constructor and Finalize
+    static_assert(sizeof(limb_t) == 4 || sizeof(limb_t) == 8, "bad size for limb_t");
+
+    void Multiply(const Num3072& a);
+    void Divide(const Num3072& a);
+    void SetToOne();
+    void Square();
+
+    Num3072() { this->SetToOne(); };
+
+    SERIALIZE_METHODS(Num3072, obj)
+    {
+        for (auto& limb : obj.limbs) {
+            READWRITE(limb);
+        }
+    }
+};
+
+/** A class representing MuHash sets
+ *
+ * MuHash is a hashing algorithm that supports adding set elements in any
+ * order but also deleting in any order. As a result, it can maintain a
+ * running sum for a set of data as a whole, and add/remove when data
+ * is added to or removed from it. A downside of MuHash is that computing
+ * an inverse is relatively expensive. This is solved by representing
+ * the running value as a fraction, and multiplying added elements into
+ * the numerator and removed elements into the denominator. Only when the
+ * final hash is desired, a single modular inverse and multiplication is
+ * needed to combine the two. The combination is also run on serialization
+ * to allow for space-efficient storage on disk.
+ *
+ * As the update operations are also associative, H(a)+H(b)+H(c)+H(d) can
+ * in fact be computed as (H(a)+H(b)) + (H(c)+H(d)). This implies that
+ * all of this is perfectly parallellizable: each thread can process an
+ * arbitrary subset of the update operations, allowing them to be
+ * efficiently combined later.
+ *
+ * Muhash does not support checking if an element is already part of the
+ * set. That is why this class does not enforce the use of a set as the
+ * data it represents because there is no efficient way to do so.
+ * It is possible to add elements more than once and also to remove
+ * elements that have not been added before. However, this implementation
+ * is intended to represent a set of elements.
+ *
+ * See also https://cseweb.ucsd.edu/~mihir/papers/inchash.pdf and
+ * https://lists.linuxfoundation.org/pipermail/bitcoin-dev/2017-May/014337.html.
+ */
+class MuHash3072
+{
+private:
+    static constexpr size_t BYTE_SIZE = 384;
+
+    Num3072 m_numerator;
+    Num3072 m_denominator;
+
+    Num3072 ToNum3072(Span<const unsigned char> in);
+
+public:
+    /* The empty set. */
+    MuHash3072() noexcept {};
+
+    /* A singleton with variable sized data in it. */
+    explicit MuHash3072(Span<const unsigned char> in) noexcept;
+
+    /* Insert a single piece of data into the set. */
+    MuHash3072& Insert(Span<const unsigned char> in) noexcept;
+
+    /* Remove a single piece of data from the set. */
+    MuHash3072& Remove(Span<const unsigned char> in) noexcept;
+
+    /* Multiply (resulting in a hash for the union of the sets) */
+    MuHash3072& operator*=(const MuHash3072& mul) noexcept;
+
+    /* Divide (resulting in a hash for the difference of the sets) */
+    MuHash3072& operator/=(const MuHash3072& div) noexcept;
+
+    /* Finalize into a 32-byte hash. Does not change this object's value. */
+    void Finalize(uint256& out) noexcept;
+
+    SERIALIZE_METHODS(MuHash3072, obj)
+    {
+        READWRITE(obj.m_numerator);
+        READWRITE(obj.m_denominator);
+    }
+};
+
+#endif // BITCOIN_CRYPTO_MUHASH_H
diff --git a/src/test/crypto_tests.cpp b/src/test/crypto_tests.cpp
index 0ad5066603..7358b246b6 100644
--- a/src/test/crypto_tests.cpp
+++ b/src/test/crypto_tests.cpp
@@ -14,7 +14,9 @@
 #include <crypto/sha256.h>
 #include <crypto/sha3.h>
 #include <crypto/sha512.h>
+#include <crypto/muhash.h>
 #include <random.h>
+#include <streams.h>
 #include <test/util/setup_common.h>
 #include <util/strencodings.h>
 
@@ -857,4 +859,92 @@ BOOST_AUTO_TEST_CASE(sha3_256_tests)
     TestSHA3_256("72c57c359e10684d0517e46653a02d18d29eff803eb009e4d5eb9e95add9ad1a4ac1f38a70296f3a369a16985ca3c957de2084cdc9bdd8994eb59b8815e0debad4ec1f001feac089820db8becdaf896aaf95721e8674e5d476b43bd2b873a7d135cd685f545b438210f9319e4dcd55986c85303c1ddf18dc746fe63a409df0a998ed376eb683e16c09e6e9018504152b3e7628ef350659fb716e058a5263a18823d2f2f6ee6a8091945a48ae1c5cb1694cf2c1fe76ef9177953afe8899cfa2b7fe0603bfa3180937dadfb66fbbdd119bbf8063338aa4a699075a3bfdbae8db7e5211d0917e9665a702fc9b0a0a901d08bea97654162d82a9f05622b060b634244779c33427eb7a29353a5f48b07cbefa72f3622ac5900bef77b71d6b314296f304c8426f451f32049b1f6af156a9dab702e8907d3cd72bb2c50493f4d593e731b285b70c803b74825b3524cda3205a8897106615260ac93c01c5ec14f5b11127783989d1824527e99e04f6a340e827b559f24db9292fcdd354838f9339a5fa1d7f6b2087f04835828b13463dd40927866f16ae33ed501ec0e6c4e63948768c5aeea3e4f6754985954bea7d61088c44430204ef491b74a64bde1358cecb2cad28ee6a3de5b752ff6a051104d88478653339457ac45ba44cbb65f54d1969d047cda746931d5e6a8b48e211416aefd5729f3d60b56b54e7f85aa2f42de3cb69419240c24e67139a11790a709edef2ac52cf35dd0a08af45926ebe9761f498ff83bfe263d6897ee97943a4b982fe3404ef0b4a45e06113c60340e0664f14799bf59cb4b3934b465fabefd87155905ee5309ba41e9e402973311831ea600b16437f71df39ee77130490c4d0227e5d1757fdc66af3ae6b9953053ed9aafca0160209858a7d4dd38fe10e0cb153672d08633ed6c54977aa0a6e67f9ff2f8c9d22dd7b21de08192960fd0e0da68d77c8d810db11dcaa61c725cd4092cbff76c8e1debd8d0361bb3f2e607911d45716f53067bdc0d89dd4889177765166a424e9fc0cb711201099dda213355e6639ac7eb86eca2ae0ab38b7f674f37ef8a6fcca1a6f52f55d9e1dcd631d2c3c82bba129172feb991d5af51afecd9d61a88b6832e4107480e392aed61a8644f551665ebff6b20953b635737a4f895e429fddcfe801f606fbda74b3bf6f5767d0fac14907fcfd0aa1d4c11b9e91b01d68052399b51a29f1ae6acd965109977c14a555cbcbd21ad8cb9f8853506d4bc21c01e62d61d7b21be1b923be54914e6b0a7ca84dd11f1159193e1184568a6134a6bbadf5b4df986edcf2019390ae841cfaa44435e28ce877d3dae4177992fa5d4e5c005876dbe3d1e63bec7dcc0942762b48b1ecc6c1a918409a8a72812a1e245c0c67be6e729c2b49bc6ee4d24a8f63e78e75db45655c26a9a78aff36fcd67117f26b8f654dca664b9f0e30681874cb749e1a692720078856286c2560b0292cc837933423147569350955c9571bf8941ba128fd339cb4268f46b94bc6ee203eb7026813706ea51c4f24c91866fc23a724bf2501327e6ae89c29f8db315dc28d2c7c719514036367e018f4835f63fdecd71f9bdced7132b6c4f8b13c69a517026fcd3622d67cb632320d5e7308f78f4b7cea11f6291b137851dc6cd6366f2785c71c3f237f81a7658b2a8d512b61e0ad5a4710b7b124151689fcb2116063fbff7e9115fed7b93de834970b838e49f8f8ba5f1f874c354078b5810a55ae289a56da563f1da6cd80a3757d6073fa55e016e45ac6cec1f69d871c92fd0ae9670c74249045e6b464787f9504128736309fed205f8df4d90e332908581298d9c75a3fa36ab0c3c9272e62de53ab290c803d67b696fd615c260a47bffad16746f18ba1a10a061bacbea9369693b3c042eec36bed289d7d12e52bca8aa1c2dff88ca7816498d25626d0f1e106ebb0b4a12138e00f3df5b1c2f49d98b1756e69b641b7c6353d99dbff050f4d76842c6cf1c2a4b062fc8e6336fa689b7c9d5c6b4ab8c15a5c20e514ff070a602d85ae52fa7810c22f8eeffd34a095b93342144f7a98d024216b3d68ed7bea047517bfcd83ec83febd1ba0e5858e2bdc1d8b1f7b0f89e90ccc432a3f930cb8209462e64556c5054c56ca2a85f16b32eb83a10459d13516faa4d23302b7607b9bd38dab2239ac9e9440c314433fdfb3ceadab4b4f87415ed6f240e017221f3b5f7ac196cdf54957bec42fe6893994b46de3d27dc7fb58ca88feb5b9e79cf20053d12530ac524337b22a3629bea52f40b06d3e2128f32060f9105847daed81d35f20e2002817434659baff64494c5b5c7f9216bfda38412a0f70511159dc73bb6bae1f8eaa0ef08d99bcb31f94f6be12c29c83df45926430b366c99fca3270c15fc4056398fdf3135b7779e3066a006961d1ac0ad1c83179ce39e87a96b722ec23aabc065badf3e188347a360772ca6a447abac7e6a44f0d4632d52926332e44a0a86bff5ce699fd063bdda3ffd4c41b53ded49fecec67f40599b934e16e3fd1bc063ad7026f8d71bfd4cbaf56599586774723194b692036f1b6bb242e2ffb9c600b5215b412764599476ce475c9e5b396fbcebd6be323dcf4d0048077400aac7500db41dc95fc7f7edbe7c9c2ec5ea89943fe13b42217eef530bbd023671509e12dfce4e1c1c82955d965e6a68aa66f6967dba48feda572db1f099d9a6dc4bc8edade852b5e824a06890dc48a6a6510ecaf8cf7620d757290e3166d431abecc624fa9ac2234d2eb783308ead45544910c633a94964b2ef5fbc409cb8835ac4147d384e12e0a5e13951f7de0ee13eafcb0ca0c04946d7804040c0a3cd088352424b097adb7aad1ca4495952f3e6c0158c02d2bcec33bfda69301434a84d9027ce02c0b9725dad118", "d894b86261436362e64241e61f6b3e6589daf64dc641f60570c4c0bf3b1f2ca3");
 }
 
+static MuHash3072 FromInt(unsigned char i) {
+    unsigned char tmp[32] = {i, 0};
+    return MuHash3072(tmp);
+}
+
+BOOST_AUTO_TEST_CASE(muhash_tests)
+{
+    uint256 out;
+
+    for (int iter = 0; iter < 10; ++iter) {
+        uint256 res;
+        int table[4];
+        for (int i = 0; i < 4; ++i) {
+            table[i] = g_insecure_rand_ctx.randbits(3);
+        }
+        for (int order = 0; order < 4; ++order) {
+            MuHash3072 acc;
+            for (int i = 0; i < 4; ++i) {
+                int t = table[i ^ order];
+                if (t & 4) {
+                    acc /= FromInt(t & 3);
+                } else {
+                    acc *= FromInt(t & 3);
+                }
+            }
+            acc.Finalize(out);
+            if (order == 0) {
+                res = out;
+            } else {
+                BOOST_CHECK(res == out);
+            }
+        }
+
+        MuHash3072 x = FromInt(g_insecure_rand_ctx.randbits(4)); // x=X
+        MuHash3072 y = FromInt(g_insecure_rand_ctx.randbits(4)); // x=X, y=Y
+        MuHash3072 z; // x=X, y=Y, z=1
+        z *= x; // x=X, y=Y, z=X
+        z *= y; // x=X, y=Y, z=X*Y
+        y *= x; // x=X, y=Y*X, z=X*Y
+        z /= y; // x=X, y=Y*X, z=1
+        z.Finalize(out);
+
+        uint256 out2;
+        MuHash3072 a;
+        a.Finalize(out2);
+
+        BOOST_CHECK_EQUAL(out, out2);
+    }
+
+    MuHash3072 acc = FromInt(0);
+    acc *= FromInt(1);
+    acc /= FromInt(2);
+    acc.Finalize(out);
+    BOOST_CHECK_EQUAL(out, uint256S("10d312b100cbd32ada024a6646e40d3482fcff103668d2625f10002a607d5863"));
+
+    MuHash3072 acc2 = FromInt(0);
+    unsigned char tmp[32] = {1, 0};
+    acc2.Insert(tmp);
+    unsigned char tmp2[32] = {2, 0};
+    acc2.Remove(tmp2);
+    acc2.Finalize(out);
+    BOOST_CHECK_EQUAL(out, uint256S("10d312b100cbd32ada024a6646e40d3482fcff103668d2625f10002a607d5863"));
+
+    // Test MuHash3072 serialization
+    MuHash3072 serchk = FromInt(1); serchk *= FromInt(2);
+    std::string ser_exp = "1fa093295ea30a6a3acdc7b3f770fa538eff537528e990e2910e40bbcfd7f6696b1256901929094694b56316de342f593303dd12ac43e06dce1be1ff8301c845beb15468fff0ef002dbf80c29f26e6452bccc91b5cb9437ad410d2a67ea847887fa3c6a6553309946880fe20db2c73fe0641adbd4e86edfee0d9f8cd0ee1230898873dc13ed8ddcaf045c80faa082774279007a2253f8922ee3ef361d378a6af3ddaf180b190ac97e556888c36b3d1fb1c85aab9ccd46e3deaeb7b7cf5db067a7e9ff86b658cf3acd6662bbcce37232daa753c48b794356c020090c831a8304416e2aa7ad633c0ddb2f11be1be316a81be7f7e472071c042cb68faef549c221ebff209273638b741aba5a81675c45a5fa92fea4ca821d7a324cb1e1a2ccd3b76c4228ec8066dad2a5df6e1bd0de45c7dd5de8070bdb46db6c554cf9aefc9b7b2bbf9f75b1864d9f95005314593905c0109b71f703d49944ae94477b51dac10a816bb6d1c700bafabc8bd86fac8df24be519a2f2836b16392e18036cb13e48c5c010000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000";
+    CDataStream ss_chk(SER_DISK, PROTOCOL_VERSION);
+    ss_chk << serchk;
+    BOOST_CHECK_EQUAL(ser_exp, HexStr(ss_chk.str()));
+
+    // Test MuHash3072 deserialization
+    MuHash3072 deserchk;
+    ss_chk >> deserchk;
+    uint256 out3;
+    serchk.Finalize(out);
+    deserchk.Finalize(out3);
+    BOOST_CHECK_EQUAL(HexStr(out), HexStr(out3));
+
+    // Test MuHash3072 overflow, meaning the internal data is larger than the modulus.
+    CDataStream ss_max(ParseHex("ffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffffff010000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000"), SER_DISK, PROTOCOL_VERSION);
+    MuHash3072 overflowchk;
+    ss_max >> overflowchk;
+
+    uint256 out4;
+    overflowchk.Finalize(out4);
+    BOOST_CHECK_EQUAL(HexStr(out4), "3a31e6903aff0de9f62f9a9f7f8b861de76ce2cda09822b90014319ae5dc2271");
+}
+
 BOOST_AUTO_TEST_SUITE_END()
diff --git a/src/test/fuzz/crypto.cpp b/src/test/fuzz/crypto.cpp
index 4783cc1c43..9668b84e7b 100644
--- a/src/test/fuzz/crypto.cpp
+++ b/src/test/fuzz/crypto.cpp
@@ -4,6 +4,7 @@
 
 #include <crypto/hmac_sha256.h>
 #include <crypto/hmac_sha512.h>
+#include <crypto/muhash.h>
 #include <crypto/ripemd160.h>
 #include <crypto/sha1.h>
 #include <crypto/sha256.h>
@@ -35,6 +36,7 @@ FUZZ_TARGET(crypto)
     CSHA512 sha512;
     SHA3_256 sha3;
     CSipHasher sip_hasher{fuzzed_data_provider.ConsumeIntegral<uint64_t>(), fuzzed_data_provider.ConsumeIntegral<uint64_t>()};
+    MuHash3072 muhash;
 
     while (fuzzed_data_provider.ConsumeBool()) {
         switch (fuzzed_data_provider.ConsumeIntegralInRange<int>(0, 2)) {
@@ -60,6 +62,12 @@ FUZZ_TARGET(crypto)
             (void)Hash(data);
             (void)Hash160(data);
             (void)sha512.Size();
+
+            if (fuzzed_data_provider.ConsumeBool()) {
+                muhash *= MuHash3072(data);
+            } else {
+                muhash /= MuHash3072(data);
+            }
             break;
         }
         case 1: {
@@ -70,10 +78,11 @@ FUZZ_TARGET(crypto)
             (void)sha256.Reset();
             (void)sha3.Reset();
             (void)sha512.Reset();
+            muhash = MuHash3072();
             break;
         }
         case 2: {
-            switch (fuzzed_data_provider.ConsumeIntegralInRange<int>(0, 9)) {
+            switch (fuzzed_data_provider.ConsumeIntegralInRange<int>(0, 10)) {
             case 0: {
                 data.resize(CHash160::OUTPUT_SIZE);
                 hash160.Finalize(data);
@@ -124,6 +133,11 @@ FUZZ_TARGET(crypto)
                 sha3.Finalize(data);
                 break;
             }
+            case 10: {
+                uint256 out;
+                muhash.Finalize(out);
+                break;
+            }
             }
             break;
         }
diff --git a/src/test/fuzz/muhash.cpp b/src/test/fuzz/muhash.cpp
new file mode 100644
index 0000000000..87359c9d39
--- /dev/null
+++ b/src/test/fuzz/muhash.cpp
@@ -0,0 +1,53 @@
+// Copyright (c) 2020 The Bitcoin Core developers
+// Distributed under the MIT software license, see the accompanying
+// file COPYING or http://www.opensource.org/licenses/mit-license.php.
+
+#include <crypto/muhash.h>
+#include <test/fuzz/FuzzedDataProvider.h>
+#include <test/fuzz/fuzz.h>
+#include <test/fuzz/util.h>
+
+#include <vector>
+
+void test_one_input(const std::vector<uint8_t>& buffer)
+{
+    FuzzedDataProvider fuzzed_data_provider{buffer.data(), buffer.size()};
+    std::vector<uint8_t> data = ConsumeRandomLengthByteVector(fuzzed_data_provider);
+    std::vector<uint8_t> data2 = ConsumeRandomLengthByteVector(fuzzed_data_provider);
+    if (data.empty()) {
+        data.resize(fuzzed_data_provider.ConsumeIntegralInRange<size_t>(1, 4096), fuzzed_data_provider.ConsumeIntegral<uint8_t>());
+    }
+    if (data2.empty()) {
+        data2.resize(fuzzed_data_provider.ConsumeIntegralInRange<size_t>(1, 4096), fuzzed_data_provider.ConsumeIntegral<uint8_t>());
+    }
+
+    data = ConsumeRandomLengthByteVector(fuzzed_data_provider);
+    data2 = ConsumeRandomLengthByteVector(fuzzed_data_provider);
+
+    MuHash3072 muhash;
+
+    // Test that MuHash result is consistent independent of order of operations
+    muhash.Insert(data);
+    muhash.Insert(data2);
+
+    uint256 out;
+    muhash.Finalize(out);
+
+    muhash = MuHash3072();
+    muhash.Insert(data2);
+    muhash.Insert(data);
+
+    uint256 out2;
+    muhash.Finalize(out2);
+
+    assert(out == out2);
+
+    // Test that removing all added elements brings the object back to it's initial state
+    muhash /= muhash;
+    muhash.Finalize(out);
+
+    MuHash3072 muhash2;
+    muhash2.Finalize(out2);
+
+    assert(out == out2);
+}