模数不超过 long long 范围时的快速乘

upd: 在伟大的 CCF 允许 int128 的使用之后本文似乎没有用武之地了，那么就把它当成时代的眼泪，留作纪念吧。

笔者的话：使用前请确保评测系统的long double严格为16B ！

模数不在 int 范围内的乘法在 OI 中运用广泛，例如Millar-Rabin，Pollard-Rho等等。这样的乘法，直接乘惨遭 long long 溢出， WA 声一片。冒险运用 __int128 ，一朝事发，倾家荡产。误食 $O(\log P)$ 的龟速乘，人傻常数大， TLE 也无处申诉。这个时候， $O(1)$ 的快速乘就能派上用场。

快速乘运用真正的long double解决这样的一个问题：给定 $x,y,P$ 满足 $0\le x,y<P$ ，计算 $xy\bmod P$ 。这里 $P< 2^{63}$ 。注意 $P$ 不一定是质数。

首先分析取模， $xy\bmod P=xy-\lfloor\frac{xy}{P}\rfloor P=(xy-\lfloor\frac{xy}{P}\rfloor P) \bmod 2^{64}$ 。所以，算出 $\lfloor\frac{xy}{P}\rfloor$ 后可以巧妙地借助自然溢出计算左右两个乘法和中间的减法，得到该式对 $2^{64}$ 取模的结果，也就得到了该式。

于是我们只剩下了一个难点就是计算 $\lfloor\frac{xy}{P}\rfloor$ 。这肯定不能先乘再除，只能先除再乘再取整。考虑搬出16B 的 long double ，先除再乘。

既然运用了 long double ，无疑这会有精度误差。极端情况是 $\frac{x}{P}$ 的第一个有效位在小数点后第一位，那么 long double 从第六十五位开始出错，因此误差范围在 $(-2^{-64},2^{-64})$ ，乘 $y$ 之后就为 $(-\frac{1}{2},\frac{1}{2})$ 。运用常见技巧将其加上0.5L，误差范围变为 $(0,1)$ 。于是，取整时误差为 $0,1$ 之一，乘 $-P$ 以后，最终误差为 $0,-P$ 之一。也就是说，设答案为 $a$ ，最后的计算结果 $r$ 为 $a-P,a$ 之一。

由于 P 在 long long 内，所以当 $0\le r<P$ 时其为第二类，直接返回 $r$ ，否则为第一类，返回 $r-P$ ，就完成了快速乘计算。时间比较充裕的话直接返回 $(r+P)\bmod P$ 即可。

经实测，在笔者笔记本电脑上，开启O2进行1e8次此数据范围内的乘法运算，快速乘耗时 6.248s ，而龟速乘耗时 55.043s ，并且计算结果相同，足以见得快速乘的优越性。

附代码：

快速乘

typedef unsigned long long ll;
ll mul(ll x,ll y,ll P){ll z=x*y-(ll)((long double)x/P*y+0.5L)*P;return z<P?z:z+P;}

 

龟速乘

typedef unsigned long long ll;
ll mul(ll x,ll y,ll P){ll z=0;for(;y;(x<<=1)>=P&&(x-=P),y>>=1)y&1&&((z+=x)>=P&&(z-=P));return z;}