O(1) 快速乘

#include<iostream>
using namespace std;
using ll = long long;
using ld = long double;

ll mul(ll a, ll b, ll p)
{
    ll c = (ld)a * b / p;
    return a * b - c * p;
}

int main()
{
    ll a, b, p; cin >> a >> b >> p;
    cout << mul(a, b, p) << '\n';
}

这是根据蓝书的快速乘修改而来。

以下是蓝书标答：

typedef unsigned long long ull;
ull mul(ull a, ull b, ull p)
{
    a %= p, b %= p;
    ull c = (long double)a * b / p;
    ull x = a * b, y = c * p;
    long long ans = (long long)(x % p) - (long long)(y % p);
    if(ans < 0) ans += p;
    return ans;
}

让咱推理一下：$$a*b\ mod\ p \Rightarrow\ a*b-\lfloor a*b/p\rfloor *p$$
这看起来不难理解。
设$c=\lfloor a*b/p\rfloor$：$$a*b\ mod\ p\Rightarrow\ a*b-c*p$$
这看起来也符合逻辑。
那：
$$a*b\ mod\ p<2^{64}\Rightarrow\ a*b-c*p<2^{64}$$
$$a*b-c*p\Rightarrow (a*b-c*p)\ mod\ 2^{64}$$
看起来就更无法击破了。
再进一步：$$(a*b-c*p)\ mod\ 2^{64}\Rightarrow (a*b\ mod\ 2^{64})-(c*p\ mod\ 2^{64})$$
好像也毫无破绽。

那这样就不会爆了吗？
那确实爆了，只不过这个叫对 $2^{64}$ 自然取模。

那$c=\lfloor a*b/p\rfloor$的时候还没推到 $2^{64}$ 啊？
那确实没到，所以用了long double来保证$a*b$不会爆。

那就有疑问了，对ull自然取模可以，难道ll就不行了吗？那自然可以。
而且我们可以推出如下等价：
$$a*b\ mod\ p\Leftrightarrow a*b-c*p\Leftrightarrow (a*b-c*p)\ mod\ 2^{63}-1$$
其中 $2^{63}-1$ 是long long的范围。
考虑到溢出可以自动取模，取模这一步也不必写了。
可以推出以下答案：

using ll = long long;
using ld = long double;
ll mul(ll a, ll b, ll p)
{
    ll c = (ld)a * b / p;
    return a * b - c * p;
}