diff --git a/src/crypto/muhash.cpp b/src/crypto/muhash.cpp
index d94d77ec2..ef3e2cbd9 100644
--- a/src/crypto/muhash.cpp
+++ b/src/crypto/muhash.cpp
@@ -1,384 +1,384 @@
 // 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;
 /**
  * 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 whether 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).
 
     // p[i] = a^(2^(2^i)-1)
     Num3072 p[12];
 
     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]);
     }
 
     Num3072 out;
 
     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::Num3072(const uint8_t (&data)[BYTE_SIZE]) {
     for (int i = 0; i < LIMBS; ++i) {
         if (sizeof(limb_t) == 4) {
             this->limbs[i] = ReadLE32(data + 4 * i);
         } else if (sizeof(limb_t) == 8) {
             this->limbs[i] = ReadLE64(data + 8 * i);
         }
     }
 }
 
 void Num3072::ToBytes(uint8_t (&out)[BYTE_SIZE]) {
     for (int i = 0; i < LIMBS; ++i) {
         if (sizeof(limb_t) == 4) {
             WriteLE32(out + i * 4, this->limbs[i]);
         } else if (sizeof(limb_t) == 8) {
             WriteLE64(out + i * 8, this->limbs[i]);
         }
     }
 }
 
 Num3072 MuHash3072::ToNum3072(Span<const uint8_t> in) {
     uint8_t tmp[Num3072::BYTE_SIZE];
 
     uint256 hashed_in = (CHashWriter(SER_DISK, 0) << in).GetSHA256();
     ChaCha20(hashed_in.data(), hashed_in.size())
         .Keystream(tmp, Num3072::BYTE_SIZE);
     Num3072 out{tmp};
 
     return out;
 }
 
 MuHash3072::MuHash3072(Span<const uint8_t> in) noexcept {
     m_numerator = ToNum3072(in);
 }
 
 void MuHash3072::Finalize(uint256 &out) noexcept {
     m_numerator.Divide(m_denominator);
     // Needed to keep the MuHash object valid
     m_denominator.SetToOne();
 
     uint8_t data[Num3072::BYTE_SIZE];
     m_numerator.ToBytes(data);
 
     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 uint8_t> in) noexcept {
     m_numerator.Multiply(ToNum3072(in));
     return *this;
 }
 
 MuHash3072 &MuHash3072::Remove(Span<const uint8_t> in) noexcept {
-    m_numerator.Divide(ToNum3072(in));
+    m_denominator.Multiply(ToNum3072(in));
     return *this;
 }