以下讨论均基于C/C++。

1. 问题引入

最近做了几道有关数学的题目，然后要用到这些较大整数的乘法（比如说NOI 2018 屠龙勇士中 1012 10 12 级别的 pi p i 相乘，还有直接上到 1018 10 18 级别的快速幂），这些在刚写代码的时候不容易想到问题，发现溢出了之后才想起这些数太大了。

这类问题虽然最终答案需要取模，但是当 a,b,p≤1018 a , b , p ≤ 10 18 时，我们没有办法直接计算 a×bmodp a × b mod p 的值—— a×b a × b 会溢出。所以我们必须想些办法，不可能为了这只大一点点的东西写一个高精度。

对于这些大数的乘法，我先讲一下网络上的普遍做法。

2. 解决方法

2.0 暴力相加法

应该大家都知道我喜欢从 1 1 开始编号，所以你看到这个方法的序号是 2.0 2.0 你就知道这个算法没有什么实际意义。

（算了我还是讲一下吧，就是对于 a×bmodp a × b mod p 重复 b b 次，每次 +a + a ，然后取模，由于一般 a a 是给到 1018 10 18 以内，所以相加不会溢出，复杂度是 O(b) O ( b ) 的，其中 b b 是 1018 10 18 级别的）

2.1 反复翻倍法

（正经讨论开始）

这个东西其实就是类似于快速幂（快速幂的一般写法通常称为反复平方法），这里只是将平方变为 ×2 × 2 。相信理解快速幂的同学一定会明白这是什么意思吧！我直接贴个代码：

inline long long qmul(long long a,long long b,long long p){
    long long res=0;
    while(b){
        if(b&1)
            res=(res+a)%p;
        a=(a<<1)%p;b>>=1;
    }
    return res;
}
//备注：这个代码是写博客时直接手打，如有错误希望在评论区帮我指出。

当然这个复杂度和快速幂的复杂度一样，是 O(log2b) O ( log 2 ⁡ b ) 的（较暴力相加法快了非常多是不是啊）。但是对于通常认为是 O(1) O ( 1 ) 的直接相乘（即C/C++的*号）还是差了不少，因此人称龟速乘，和它的祖先快速幂形成了鲜明的对比。

2.2 位运算技巧法

所以为了避免出现 n=106 n = 10 6 的大数据而导致的 O(nlog2ai) O ( n log 2 ⁡ a i ) 的复杂度被卡，
我们还是希望找到一种理论复杂度 O(1) O ( 1 ) 的方法。这是存在的。

Notes：以下讨论均默认 0≤a,b<p 0 ≤ a , b < p 。如果出现负数或是超过模数的情况，请自行预先处理。

首先，我们先明确模运算的定义：

a×bmodp=a×b−⌊a×bp⌋×p a × b mod p = a × b − ⌊ a × b p ⌋ × p

所以，当 a,b<p a , b < p 时，一定满足 ⌊a×bp⌋<p ⌊ a × b p ⌋ < p 。

然后我们可以用long double来存储 a×bp a × b p 的值（浮点数形式），根据浮点数的运算规则，当精度不够用的时候，舍弃末尾也就是小数点后的几位，因此使用long double即可满足要求（现在绝大多数平台上long double给了 12 12 或 16 16 个字节，因此有效位数是 18−19 18 − 19 位，可以使用这种方法，对于那些老旧的机器例如NOIP评测机在不确定的情况下还是使用龟速乘吧）。

接下来我们的得到的是一个拖着好几位小数的浮点数，然后我们把它扔进一个long long里（不妨记为 c c ，那么 c=⌊a×bp⌋ c = ⌊ a × b p ⌋ ），以达到下取整的目的。（严重警告：请不要直接使用long double将原式算出来，否则可能产生严重误差或发生溢出。请只在这一步使用）。然后我们就可以把答案看成 a×b−c×p a × b − c × p 了。

慢着，这里long long类型的 a a 和 b b 相乘不会溢出吗？

是的， a×b a × b 会溢出， c×p c × p 也会溢出，但是无论如何，他们之间的差值总是在 [0,p) [ 0 , p ) 内，所以我们只需要关心最后的十几位即可，那些溢出的只是舍去了高位，对低位没有影响，因此答案是正确的。

综上，我们只是根据模运算的定义，然后结合了long double字节长但是不精确的特点，避开了 64 64 位的限制算出答案。技巧性非常强，但值得一学。

下面给出代码：

inline long long qmul(long long a,long long b,long long p){
    a=(a%p+p)%p;b=(b%p+p)%p;//顺带处理了负数和过大的情况
    long long c=a*(long double)b/p;//这里的c是准确的下取整结果
    long long ans=a*b-c*p;
    if(ans<0)//如果不在[0,p)之间则调整一下
        ans+=p;
    else if(ans>=p)
        ans-=p;
    return ans;
}

时间复杂度 O(1) O ( 1 ) 。

3. 总结

技巧性很强，但要注意细节（而且小数据难以查错）。