gistfile1.cpp

// Computes cR such that cR = aR * bR (mod n) with R = 2**64.
inline uint64_t montmul64(uint64_t aR, uint64_t bR, uint64_t n, uint64_t nneginv) {
    uint64_t Th, Tl, m, mnh, mnl, th;

    std::tie(Th, Tl) = op::mulu64(aR, bR);
    m = Tl * nneginv;
    std::tie(mnh, mnl) = op::mulu64(m, n);

    bool lc = Tl + mnl < Tl;
    th = Th + mnh + lc;
    bool hc = (th < Th) || (th == Th && lc);

    if (hc > 0 || th >= n) th = th - n;

    return th;
}