笔者的话：使用前请确保评测系统的long double严格为16B ！
模数不在 int 范围内的乘法在 OI 中运用广泛，例如Millar-Rabin，Pollard-Rho等等。这样的乘法，直接乘惨遭 long long 溢出， WA 声一片。冒险运用 __int128 ，一朝事发，倾家荡产。误食 $O(log P)$ 的龟速乘，人傻常数大， TLE 也无处申诉。这个时候， $O(1)$ 的快速乘就能派上用场。
快速乘运用真正的long double解决这样的一个问题：给定 $x,y,P$ 满足 $0le x,y<P$ ，计算 $xymod P$ 。这里 $P< 2^{63}$ 。注意 $P$ 不一定是质数。
首先分析取模， $xymod P=xy-lfloorfrac{xy}{P} floor P=(xy-lfloorfrac{xy}{P} floor P) mod 2^{64}$ 。所以，算出 $lfloorfrac{xy}{P} floor$ 后可以巧妙地借助自然溢出计算左右两个乘法和中间的减法，得到该式对 $2^{64}$ 取模的结果，也就得到了该式。
于是我们只剩下了一个难点就是计算 $lfloorfrac{xy}{P} floor$ 。这肯定不能先乘再除，只能先除再乘再取整。考虑搬出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 内，所以当 $0le r<P$ 时其为第二类，直接返回 $r$ ，否则为第一类，返回 $r-P$ ，就完成了快速乘计算。时间比较充裕的话直接返回 $(r+P)mod 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;}