muhash.cpp

Go to the documentation of this file.

// 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 <util/check.h>
 
#include <cassert>
#include <cstdio>
#include <limits>
 
namespace {
 
using limb_t = Num3072::limb_t;
using signed_limb_t = Num3072::signed_limb_t;
using double_limb_t = Num3072::double_limb_t;
using signed_double_limb_t = Num3072::signed_double_limb_t;
constexpr int LIMB_SIZE = Num3072::LIMB_SIZE;
constexpr int SIGNED_LIMB_SIZE = Num3072::SIGNED_LIMB_SIZE;
constexpr int LIMBS = Num3072::LIMBS;
constexpr int SIGNED_LIMBS = Num3072::SIGNED_LIMBS;
constexpr int FINAL_LIMB_POSITION = 3072 / SIGNED_LIMB_SIZE;
constexpr int FINAL_LIMB_MODULUS_BITS = 3072 % SIGNED_LIMB_SIZE;
constexpr limb_t MAX_LIMB = (limb_t)(-1);
constexpr limb_t MAX_SIGNED_LIMB = (((limb_t)1) << SIGNED_LIMB_SIZE) - 1;
constexpr limb_t MAX_PRIME_DIFF = 1103717;
constexpr limb_t MODULUS_INVERSE = limb_t(0x70a1421da087d93);
 
inline void extract3(limb_t &c0, limb_t &c1, limb_t &c2, limb_t &n) {
    n = c0;
    c0 = c1;
    c1 = c2;
    c2 = 0;
}
 
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;
}
 
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;
}
 
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;
}
 
} // namespace
 
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]);
    }
}
 
namespace {
struct Num3072Signed {
    signed_limb_t limbs[SIGNED_LIMBS];
 
    Num3072Signed() { memset(limbs, 0, sizeof(limbs)); }
 
    void FromNum3072(const Num3072 &in) {
        double_limb_t c = 0;
        int b = 0, outpos = 0;
        for (int i = 0; i < LIMBS; ++i) {
            c += double_limb_t{in.limbs[i]} << b;
            b += LIMB_SIZE;
            while (b >= SIGNED_LIMB_SIZE) {
                limbs[outpos++] = limb_t(c) & MAX_SIGNED_LIMB;
                c >>= SIGNED_LIMB_SIZE;
                b -= SIGNED_LIMB_SIZE;
            }
        }
        Assume(outpos == SIGNED_LIMBS - 1);
        limbs[SIGNED_LIMBS - 1] = c;
        c >>= SIGNED_LIMB_SIZE;
        Assume(c == 0);
    }
 
    void ToNum3072(Num3072 &out) const {
        double_limb_t c = 0;
        int b = 0, outpos = 0;
        for (int i = 0; i < SIGNED_LIMBS; ++i) {
            c += double_limb_t(limbs[i]) << b;
            b += SIGNED_LIMB_SIZE;
            if (b >= LIMB_SIZE) {
                out.limbs[outpos++] = c;
                c >>= LIMB_SIZE;
                b -= LIMB_SIZE;
            }
        }
        Assume(outpos == LIMBS);
        Assume(c == 0);
    }
 
    void Normalize(bool negate) {
        // Add modulus if this was negative. This brings the range of *this to
        // 1-2^3072..2^3072-1. -1 if this is negative; 0 otherwise
        signed_limb_t cond_add = limbs[SIGNED_LIMBS - 1] >> (LIMB_SIZE - 1);
        limbs[0] += signed_limb_t(-MAX_PRIME_DIFF) & cond_add;
        limbs[FINAL_LIMB_POSITION] +=
            (signed_limb_t(1) << FINAL_LIMB_MODULUS_BITS) & cond_add;
        // Next negate all limbs if negate was set. This does not change the
        // range of *this. -1 if this negate is true; 0 otherwise
        signed_limb_t cond_negate = -signed_limb_t(negate);
        for (int i = 0; i < SIGNED_LIMBS; ++i) {
            limbs[i] = (limbs[i] ^ cond_negate) - cond_negate;
        }
        // Perform carry (make all limbs except the top one be in range
        // 0..2^SIGNED_LIMB_SIZE-1).
        for (int i = 0; i < SIGNED_LIMBS - 1; ++i) {
            limbs[i + 1] += limbs[i] >> SIGNED_LIMB_SIZE;
            limbs[i] &= MAX_SIGNED_LIMB;
        }
        // Again add modulus if *this was negative. This brings the range of
        // *this to 0..2^3072-1. -1 if this is negative; 0 otherwise
        cond_add = limbs[SIGNED_LIMBS - 1] >> (LIMB_SIZE - 1);
        limbs[0] += signed_limb_t(-MAX_PRIME_DIFF) & cond_add;
        limbs[FINAL_LIMB_POSITION] +=
            (signed_limb_t(1) << FINAL_LIMB_MODULUS_BITS) & cond_add;
        // Perform another carry. Now all limbs are in range
        // 0..2^SIGNED_LIMB_SIZE-1.
        for (int i = 0; i < SIGNED_LIMBS - 1; ++i) {
            limbs[i + 1] += limbs[i] >> SIGNED_LIMB_SIZE;
            limbs[i] &= MAX_SIGNED_LIMB;
        }
    }
};
 
struct SignedMatrix {
    signed_limb_t u, v, q, r;
};
 
inline int CountTrailingZeroes(limb_t v) {
#ifdef HAVE_DECL___BUILTIN_CTZ
    /* Use __builtin_ctz if available, and sufficient for the size of limb_t. */
    if constexpr (std::numeric_limits<limb_t>::max() <=
                  std::numeric_limits<unsigned>::max()) {
        return __builtin_ctz(v);
    }
#endif
#ifdef HAVE_DECL___BUILTIN_CTZL
    /* Use __builtin_ctzl if available, and sufficient for the size of limb_t.
     */
    if constexpr (std::numeric_limits<limb_t>::max() <=
                  std::numeric_limits<unsigned long>::max()) {
        return __builtin_ctzl(v);
    }
#endif
#ifdef HAVE_DECL___BUILTIN_CTZLL
    /* Use __builtin_ctzll if available, and sufficient for the size of limb_t.
     */
    if constexpr (std::numeric_limits<limb_t>::max() <=
                  std::numeric_limits<unsigned long long>::max()) {
        return __builtin_ctzll(v);
    }
#endif
    /* Otherwise fall back to a software implementation. */
    if constexpr (LIMB_SIZE <= 32) {
        static constexpr uint8_t debruijn[32] = {
            0,  1,  2,  24, 3,  19, 6,  25, 26, 4,  20, 10, 16, 7,  12, 26,
            31, 23, 18, 5,  21, 9,  15, 11, 30, 17, 8,  14, 29, 13, 28, 27};
        return debruijn[((v & -v) * 0x04D7651F) >> 27];
    } else {
        static_assert(LIMB_SIZE <= 64);
        static constexpr uint8_t debruijn[64] = {
            0,  1,  2,  53, 3,  7,  54, 27, 4,  38, 41, 8,  34, 55, 48, 28,
            62, 5,  39, 46, 44, 42, 22, 9,  24, 35, 59, 56, 49, 18, 29, 11,
            63, 52, 6,  26, 37, 40, 33, 47, 61, 45, 43, 21, 23, 58, 17, 10,
            51, 25, 36, 32, 60, 20, 57, 16, 50, 31, 19, 15, 30, 14, 13, 12};
        return debruijn[((v & -v) * 0x022FDD63CC95386D) >> 58];
    }
}
 
inline limb_t ComputeDivstepMatrix(signed_limb_t eta, limb_t f, limb_t g,
                                   SignedMatrix &out) {
    static const uint8_t NEGINV256[128] = {
        0xFF, 0x55, 0x33, 0x49, 0xC7, 0x5D, 0x3B, 0x11, 0x0F, 0xE5, 0xC3, 0x59,
        0xD7, 0xED, 0xCB, 0x21, 0x1F, 0x75, 0x53, 0x69, 0xE7, 0x7D, 0x5B, 0x31,
        0x2F, 0x05, 0xE3, 0x79, 0xF7, 0x0D, 0xEB, 0x41, 0x3F, 0x95, 0x73, 0x89,
        0x07, 0x9D, 0x7B, 0x51, 0x4F, 0x25, 0x03, 0x99, 0x17, 0x2D, 0x0B, 0x61,
        0x5F, 0xB5, 0x93, 0xA9, 0x27, 0xBD, 0x9B, 0x71, 0x6F, 0x45, 0x23, 0xB9,
        0x37, 0x4D, 0x2B, 0x81, 0x7F, 0xD5, 0xB3, 0xC9, 0x47, 0xDD, 0xBB, 0x91,
        0x8F, 0x65, 0x43, 0xD9, 0x57, 0x6D, 0x4B, 0xA1, 0x9F, 0xF5, 0xD3, 0xE9,
        0x67, 0xFD, 0xDB, 0xB1, 0xAF, 0x85, 0x63, 0xF9, 0x77, 0x8D, 0x6B, 0xC1,
        0xBF, 0x15, 0xF3, 0x09, 0x87, 0x1D, 0xFB, 0xD1, 0xCF, 0xA5, 0x83, 0x19,
        0x97, 0xAD, 0x8B, 0xE1, 0xDF, 0x35, 0x13, 0x29, 0xA7, 0x3D, 0x1B, 0xF1,
        0xEF, 0xC5, 0xA3, 0x39, 0xB7, 0xCD, 0xAB, 0x01};
    // Coefficients of returned SignedMatrix; starts off as identity matrix. */
    limb_t u = 1, v = 0, q = 0, r = 1;
    // The number of divsteps still left.
    int i = SIGNED_LIMB_SIZE;
    while (true) {
        /* Use a sentinel bit to count zeros only up to i. */
        int zeros = CountTrailingZeroes(g | (MAX_LIMB << i));
        /* Perform zeros divsteps at once; they all just divide g by two. */
        g >>= zeros;
        u <<= zeros;
        v <<= zeros;
        eta -= zeros;
        i -= zeros;
        /* We're done once we've performed SIGNED_LIMB_SIZE divsteps. */
        if (i == 0) {
            break;
        }
        /* If eta is negative, negate it and replace f,g with g,-f. */
        if (eta < 0) {
            limb_t tmp;
            eta = -eta;
            tmp = f;
            f = g;
            g = -tmp;
            tmp = u;
            u = q;
            q = -tmp;
            tmp = v;
            v = r;
            r = -tmp;
        }
        /* eta is now >= 0. In what follows we're going to cancel out the bottom
         * bits of g. No more than i can be cancelled out (as we'd be done
         * before that point), and no more than eta+1 can be done as its sign
         * will flip once that happens.
         */
        int limit = ((int)eta + 1) > i ? i : ((int)eta + 1);
        /*
         * m is a mask for the bottom min(limit, 8) bits (our table only
         * supports 8 bits).
         */
        limb_t m = (MAX_LIMB >> (LIMB_SIZE - limit)) & 255U;
        /*
         * Find what multiple of f must be added to g to cancel its bottom
         * min(limit, 8) bits.
         */
        limb_t w = (g * NEGINV256[(f >> 1) & 127]) & m;
        /* Do so. */
        g += f * w;
        q += u * w;
        r += v * w;
    }
    out.u = (signed_limb_t)u;
    out.v = (signed_limb_t)v;
    out.q = (signed_limb_t)q;
    out.r = (signed_limb_t)r;
    return eta;
}
 
inline void UpdateDE(Num3072Signed &d, Num3072Signed &e,
                     const SignedMatrix &t) {
    const signed_limb_t u = t.u, v = t.v, q = t.q, r = t.r;
 
    /*
     * [md,me] start as zero; plus [u,q] if d is negative; plus [v,r] if e is
     * negative.
     */
    signed_limb_t sd = d.limbs[SIGNED_LIMBS - 1] >> (LIMB_SIZE - 1);
    signed_limb_t se = e.limbs[SIGNED_LIMBS - 1] >> (LIMB_SIZE - 1);
    signed_limb_t md = (u & sd) + (v & se);
    signed_limb_t me = (q & sd) + (r & se);
    /* Begin computing t*[d,e]. */
    signed_limb_t di = d.limbs[0], ei = e.limbs[0];
    signed_double_limb_t cd =
        (signed_double_limb_t)u * di + (signed_double_limb_t)v * ei;
    signed_double_limb_t ce =
        (signed_double_limb_t)q * di + (signed_double_limb_t)r * ei;
    /*
     * Correct md,me so that t*[d,e]+modulus*[md,me] has SIGNED_LIMB_SIZE zero
     * bottom bits.
     */
    md -= (MODULUS_INVERSE * limb_t(cd) + md) & MAX_SIGNED_LIMB;
    me -= (MODULUS_INVERSE * limb_t(ce) + me) & MAX_SIGNED_LIMB;
    /*
     * Update the beginning of computation for t*[d,e]+modulus*[md,me] now md,me
     * are known.
     */
    cd -= (signed_double_limb_t)1103717 * md;
    ce -= (signed_double_limb_t)1103717 * me;
    /*
     * Verify that the low SIGNED_LIMB_SIZE bits of the computation are indeed
     * zero, and then throw them away.
     */
    Assume((cd & MAX_SIGNED_LIMB) == 0);
    Assume((ce & MAX_SIGNED_LIMB) == 0);
    cd >>= SIGNED_LIMB_SIZE;
    ce >>= SIGNED_LIMB_SIZE;
    /*
     * Now iteratively compute limb i=1..SIGNED_LIMBS-2 of
     * t*[d,e]+modulus*[md,me], and store them in output limb i-1 (shifting down
     * by SIGNED_LIMB_SIZE bits). The corresponding limbs in modulus are all
     * zero, so modulus/md/me are not actually involved here.
     */
    for (int i = 1; i < SIGNED_LIMBS - 1; ++i) {
        di = d.limbs[i];
        ei = e.limbs[i];
        cd += (signed_double_limb_t)u * di + (signed_double_limb_t)v * ei;
        ce += (signed_double_limb_t)q * di + (signed_double_limb_t)r * ei;
        d.limbs[i - 1] = (signed_limb_t)cd & MAX_SIGNED_LIMB;
        cd >>= SIGNED_LIMB_SIZE;
        e.limbs[i - 1] = (signed_limb_t)ce & MAX_SIGNED_LIMB;
        ce >>= SIGNED_LIMB_SIZE;
    }
    /*
     * Compute limb SIGNED_LIMBS-1 of t*[d,e]+modulus*[md,me], and store it in
     * output limb SIGNED_LIMBS-2.
     */
    di = d.limbs[SIGNED_LIMBS - 1];
    ei = e.limbs[SIGNED_LIMBS - 1];
    cd += (signed_double_limb_t)u * di + (signed_double_limb_t)v * ei;
    ce += (signed_double_limb_t)q * di + (signed_double_limb_t)r * ei;
    cd += (signed_double_limb_t)md << FINAL_LIMB_MODULUS_BITS;
    ce += (signed_double_limb_t)me << FINAL_LIMB_MODULUS_BITS;
    d.limbs[SIGNED_LIMBS - 2] = (signed_limb_t)cd & MAX_SIGNED_LIMB;
    cd >>= SIGNED_LIMB_SIZE;
    e.limbs[SIGNED_LIMBS - 2] = (signed_limb_t)ce & MAX_SIGNED_LIMB;
    ce >>= SIGNED_LIMB_SIZE;
    /* What remains goes into output limb SINGED_LIMBS-1 */
    d.limbs[SIGNED_LIMBS - 1] = (signed_limb_t)cd;
    e.limbs[SIGNED_LIMBS - 1] = (signed_limb_t)ce;
}
 
inline void UpdateFG(Num3072Signed &f, Num3072Signed &g, const SignedMatrix &t,
                     int len) {
    const signed_limb_t u = t.u, v = t.v, q = t.q, r = t.r;
 
    signed_limb_t fi, gi;
    signed_double_limb_t cf, cg;
    /* Start computing t*[f,g]. */
    fi = f.limbs[0];
    gi = g.limbs[0];
    cf = (signed_double_limb_t)u * fi + (signed_double_limb_t)v * gi;
    cg = (signed_double_limb_t)q * fi + (signed_double_limb_t)r * gi;
    /*
     * Verify that the bottom SIGNED_LIMB_BITS bits of the result are zero, and
     * then throw them away.
     */
    Assume((cf & MAX_SIGNED_LIMB) == 0);
    Assume((cg & MAX_SIGNED_LIMB) == 0);
    cf >>= SIGNED_LIMB_SIZE;
    cg >>= SIGNED_LIMB_SIZE;
    /*
     * Now iteratively compute limb i=1..SIGNED_LIMBS-1 of t*[f,g], and store
     * them in output limb i-1 (shifting down by SIGNED_LIMB_BITS bits).
     */
    for (int i = 1; i < len; ++i) {
        fi = f.limbs[i];
        gi = g.limbs[i];
        cf += (signed_double_limb_t)u * fi + (signed_double_limb_t)v * gi;
        cg += (signed_double_limb_t)q * fi + (signed_double_limb_t)r * gi;
        f.limbs[i - 1] = (signed_limb_t)cf & MAX_SIGNED_LIMB;
        cf >>= SIGNED_LIMB_SIZE;
        g.limbs[i - 1] = (signed_limb_t)cg & MAX_SIGNED_LIMB;
        cg >>= SIGNED_LIMB_SIZE;
    }
    /*
     * What remains is limb SIGNED_LIMBS of t*[f,g]; store it as output limb
     * SIGNED_LIMBS-1.
     */
    f.limbs[len - 1] = (signed_limb_t)cf;
    g.limbs[len - 1] = (signed_limb_t)cg;
}
} // namespace
 
Num3072 Num3072::GetInverse() const {
    // Compute a modular inverse based on a variant of the safegcd algorithm:
    // - Paper: https://gcd.cr.yp.to/papers.html
    // - Inspired by this code in libsecp256k1:
    //   https://github.com/bitcoin-core/secp256k1/blob/master/src/modinv32_impl.h
    // - Explanation of the algorithm:
    //   https://github.com/bitcoin-core/secp256k1/blob/master/doc/safegcd_implementation.md
 
    // Local variables d, e, f, g:
    // - f and g are the variables whose gcd we compute (despite knowing the
    // answer is 1):
    //   - f is always odd, and initialized as modulus
    //   - g is initialized as *this (called x in what follows)
    // - d and e are the numbers for which at every step it is the case that:
    //   - f = d * x mod modulus; d is initialized as 0
    //   - g = e * x mod modulus; e is initialized as 1
    Num3072Signed d, e, f, g;
    e.limbs[0] = 1;
    // F is initialized as modulus, which in signed limb representation can be
    // expressed simply as 2^3072 + -MAX_PRIME_DIFF.
    f.limbs[0] = -MAX_PRIME_DIFF;
    f.limbs[FINAL_LIMB_POSITION] = ((limb_t)1) << FINAL_LIMB_MODULUS_BITS;
    g.FromNum3072(*this);
    int len = SIGNED_LIMBS; 
    signed_limb_t eta = -1; 
    // Perform divsteps on [f,g] until g=0 is reached, keeping (d,e)
    // synchronized with them.
    while (true) {
        // Compute transformation matrix t that represents the next
        // SIGNED_LIMB_SIZE divsteps to apply. This can be computed from just
        // the bottom limb of f and g, and eta.
        SignedMatrix t;
        eta = ComputeDivstepMatrix(eta, f.limbs[0], g.limbs[0], t);
        // Apply that transformation matrix to the full [f,g] vector.
        UpdateFG(f, g, t, len);
        // Apply that transformation matrix to the full [d,e] vector (mod
        // modulus).
        UpdateDE(d, e, t);
 
        // Check if g is zero.
        if (g.limbs[0] == 0) {
            signed_limb_t cond = 0;
            for (int j = 1; j < len; ++j) {
                cond |= g.limbs[j];
            }
            // If so, we're done.
            if (cond == 0) {
                break;
            }
        }
 
        // Check if the top limbs of both f and g are both 0 or -1.
        signed_limb_t fn = f.limbs[len - 1], gn = g.limbs[len - 1];
        signed_limb_t cond = ((signed_limb_t)len - 2) >> (LIMB_SIZE - 1);
        cond |= fn ^ (fn >> (LIMB_SIZE - 1));
        cond |= gn ^ (gn >> (LIMB_SIZE - 1));
        if (cond == 0) {
            // If so, drop the top limb, shrinking the size of f and g, by
            // propagating the sign to the previous limb.
            f.limbs[len - 2] |= (limb_t)f.limbs[len - 1] << SIGNED_LIMB_SIZE;
            g.limbs[len - 2] |= (limb_t)g.limbs[len - 1] << SIGNED_LIMB_SIZE;
            --len;
        }
    }
    // At some point, [f,g] will have been rewritten into [f',0], such that
    // gcd(f,g) = gcd(f',0). This is proven in the paper. As f started out being
    // modulus, a prime number, we know that gcd is 1, and thus f' is 1 or -1.
    Assume((f.limbs[0] & MAX_SIGNED_LIMB) == 1 ||
           (f.limbs[0] & MAX_SIGNED_LIMB) == MAX_SIGNED_LIMB);
    // As we've maintained the invariant that f = d * x mod modulus, we get d/f
    // mod modulus is the modular inverse of x we're looking for. As f is 1 or
    // -1, it is also true that d/f = d*f. Normalize d to prepare it for output,
    // while negating it if f is negative.
    d.Normalize(f.limbs[len - 1] >> (LIMB_SIZE - 1));
    Num3072 ret;
    d.ToNum3072(ret);
    return ret;
}
 
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);
 
    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{(HashWriter{} << in).GetSHA256()};
    static_assert(sizeof(tmp) % ChaCha20Aligned::BLOCKLEN == 0);
    ChaCha20Aligned{MakeByteSpan(hashed_in)}.Keystream(
        MakeWritableByteSpan(tmp));
    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 = (HashWriter{} << 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_denominator.Multiply(ToNum3072(in));
    return *this;
}