// Returns (ab)R mod n given aR mod n, bR mod n, n and -n^(-1) mod R, with R = 2^32
inline uint64_t montmul32(uint64_t a, uint64_t b, uint64_t n, uint64_t nneginv) {
uint64_t T = a*b;
uint32_t m = T*nneginv; // m = T*-n^(-1) (mod 2^32)
uint64_t t = (T + m*n) >> 32;
return t >= n ? t - n : t;
}