Barrett 取模优化
前言
今天模拟赛有人用双 \(\log\) 做法加 Barrett 优化卡过正解单 \(\log\) 都要卡常的题了，感觉很有意思，学点新东西,记录一下。
另外不知道这个算法怎么分类（因为沾点数学但又是只有写代码时才有用的）索性放到杂项里了。
Barrett 是什么
Barrett 是用来优化对同一个数的取模的，在固定模数 \(P\) 的情况下，Barrett 的速度比 C++ 自带的 % 要快得多。
Barrett 的原理是什么
首先考虑取模写成四则运算是什么，容易发现 \(a\bmod P=a-P\left\lfloor\frac{a}{P}\right\rfloor\)。也就是说，我们只要优化了算 \(\left\lfloor\frac{a}{P}\right\rfloor\) 的效率，就能加速取模。
由于 C++ 默认的除法效率并不高，考虑用估算的方法，使得 \(\frac{1}{P}\approx\frac{m}{2^k}\)，因为计算机算乘法是很快的，而 \(\frac{1}{2^k}\) 可以用右移运算解决。那么此时 \(m\approx\frac{2^k}{P}\)，不妨设 \(m=\left\lfloor\frac{2^k}{P}\right\rfloor\)，去掉下取整就是 \(m=\frac{2^k}{P}-\lambda\)，其中 \(0\le \lambda<1\)。
现在用这个近似值估算：
\[\left\lfloor\dfrac{a}{P}\right\rfloor\approx\left\lfloor\dfrac{am}{2^k}\right\rfloor=\left\lfloor a\left(\dfrac{\frac{2^k}{P}-\lambda}{2^k}\right)\right\rfloor=\left\lfloor\dfrac{a}{P}-\dfrac{a\lambda}{2^k}\right\rfloor\]
假如 \(k\) 足够大，使得 \(\dfrac{a\lambda}{2^k}<1\)，那么分类讨论一下就能得到以下式子：
\[
\left\lfloor\dfrac{a}{P}-\dfrac{a\lambda}{2^k}\right\rfloor=\begin{cases}
\left\lfloor\dfrac{a}{P}\right\rfloor, &\dfrac{a}{P}-\dfrac{a\lambda}{2^k}\ge\left\lfloor\dfrac{a}{P}\right\rfloor\\
\left\lfloor\dfrac{a}{P}\right\rfloor-1,&\dfrac{a}{P}-\dfrac{a\lambda}{2^k}<\left\lfloor\dfrac{a}{P}\right\rfloor
\end{cases}
\]
然后塞进原式里：
\[
a-P\left\lfloor\dfrac{a}{P}-\dfrac{a\lambda}{2^k}\right\rfloor=\begin{cases}
a-P\left\lfloor\dfrac{a}{P}\right\rfloor=a\bmod P, &\dfrac{a}{P}-\dfrac{a\lambda}{2^k}\ge\left\lfloor\dfrac{a}{P}\right\rfloor\\
a-P\left(\left\lfloor\dfrac{a}{P}\right\rfloor-1\right)=a\bmod P+P,&\dfrac{a}{P}-\dfrac{a\lambda}{2^k}<\left\lfloor\dfrac{a}{P}\right\rfloor
\end{cases}
\]
那么实现的时候随便特判一下即可。
现在想一下如何保证 \(\dfrac{a\lambda}{2^k}<1\)，容易发现 \(a\) 最大是 \(P-1\)，那么 \(2^k\ge P\) 即可，为了方便我们一般取 \(k=64\)（unsigned long long 的最大值加一），注意会爆 unsigned long long，所以中间要转化成 128 位无符号整数 __uint128_t。
可以用一个结构体封装一下，然后重载一下 % 符号。
参考代码
#define ull unsigned long long
#define ui128 __uint128_t
struct Barrett{
ull d;ui128 m;
void init(ull _d){
d=_d,m=(((ui128)(1)<<64)/d);
}
}mod;
ull operator %(ull a,Barrett mod){
ull w=(mod.m*a)>>64;w=a-w*mod.d;
if(w>=mod.d)w-=mod.d;return w;
}