取模优化
常见的做法是加法取模转减法，或者积累数据之后一次性取模，不过这些没有真正的加快取模运算，只是减少了取模的运算次数。
真正想要把取模做到和加减一样快（\(O(1)\)）已经被证明是不行的了，不过我们可以用 \(O(1)\) 次整除和取模，把对于常数取模优化到 \(O(1)\)。（这里认为乘法是 \(O(1)\) 的）
蒙哥马利算法
思路
上面说了，加减取模可以当做加减，也很容易实现：
inline unsigned add(unsigned x,unsigned y) { x+=y; return x>=mod?x-mod:x; }
inline unsigned sub(unsigned x,unsigned y) { return x>y?x-y:mod-(y-x); }
现在的问题就是我们计算乘法时可能会溢出很多位，不能再用减法了。蒙哥马利模乘可以把对常数的取模转化为对 \(2^k\) 的取模，如果模数是 \(2^k\)，那就可以把取模和整除转化为按位与和位移了，这显然是 \(O(1)\) 的。
假设我们的模数是 \(M\)，并且 \(M\perp 2\)。（如果不是这样，我们就不能用蒙哥马利模乘了）
再假设 \(2^{k}>M\)，显然也有 \(M\perp 2^k\)，这里的 \(2^k\) 就被叫做蒙哥马利模数。后文中，我们记 \(R=2^k\)。
显然的想法是我们进行数域转换，即 \(\mathbb{Z}_{M}\rightarrow\mathbb{Z}_{R}\rightarrow\mathbb{Z}_M\)。
首先，我们找到 \(a^\prime\equiv aR\pmod M,b^\prime\equiv bR\pmod M\)，那么就有 \(abR\equiv a^\prime b^\prime R^{-1}\pmod M\)，至此只需要使用某种方法快速算出 \(xR^{-1}\pmod M\) 即可。（我们假设了 \(a,b\in \mathbb{Z}_M\)，如果不是这样，我们也可以先用蒙哥马利算法算出 \(a\pmod M\) 和 \(b\pmod M\) 之后再计算 \(ab\pmod M\) ）
完整算法
求出 \(a^\prime\equiv aR\pmod M\)：
我们考虑预处理出 \(r=R^2\pmod M\)，那么 \(a^\prime\equiv aR\equiv arR^{-1}\pmod M\)，只需求出 \(arR^{-1}\pmod M\)，这和我们下一步要处理的问题是一样的。
求出 \(a^\prime b^\prime R^{-1}\pmod M\)：
我们原本的值域是 \(M\)，那么 \(a^\prime b^\prime\) 的值域就是 \(M^2\)，在实际运用中我们常常可以用 unsigned long long 或者 __uint128_t 之类的存下，最多写个小高精。那么这就是下一步计算 \(cR^{-1}\pmod M\) 可以完成的事了。
求出 \(cR^{-1}\pmod M\)
有 \(cR^{-1}=\dfrac{c}{R}\)，如果 \(c\) 是 \(R\) 的整数倍，那么我们可以直接位移了，否则我们可以求出 \(R\mid(c+kM)\)，那么 \(cR^{-1}\equiv\frac{c+kM}{R}\pmod M\)。
现在的问题是要求 \(k\)，稍微变换一下可以得到 \(k\equiv c(-M)^{-1}\pmod R\)，只需要预处理出 \((-M)^{-1}\pmod R\) 即可。
那么如果 \(c\in [0,MR)\)，直接位移就做完了。（代码里倒着实现了三个功能）
typedef unsigned long long ull;
unsigned R=1u<<len;
unsigned r=1ull*R*R%M,m_=pow(R-M,R>>1);
inline unsigned getK(ull c) { return c*m_>>len; }
inline unsigned MontR(ull c) { return (c+getK(c)*M)>>len; }
inline unsigned MontMul(unsigned x,unsigned y) { return MontR(1ull*x*y); }
inline unsigned MontTurn(unsigned x) { return MontMul(x,r); }
inline unsigned modMul(unsigned x,unsigned y) { return MontR(MontMul(MontTurn(x),MontTurn(y))); }
Barrett reduction
思路
运用自然溢出和位运算替代对常数整除，然后对常数取模用对常数整除实现。
比起蒙哥马利算法，Barrett reduction 对于任意常数都适用。
完整算法
预处理出 \(m,l\)，使得 \(\lfloor\dfrac nM\rfloor=\lfloor\dfrac{nm}{2^{N+l}}\rfloor\)，其中 \(2^N\) 是存储变量的类型的值域。（比如 unsigned 的 \(N\) 就取 \(32\)）
先直接给出定理：若 \(m,M,l\) 满足 \(2^N\le mM\le 2^N+2^l\)，则 \(\forall n\in [0,2^N)\)，\(\lfloor\dfrac nM\rfloor=\lfloor\dfrac{nm}{2^{N+l}}\rfloor\)。
设 \(n=qM+r,0\le r<M\) 和 \(k=mM-2^{N+l}\)，显然有 \(k\in [0,2^l]\)。
考虑
\[\begin{align*} \frac{nm}{2^{N+l}}-q&=\frac{mM}{M}\frac{n}{2^{N+l}}-q\\ &=\frac{k+2^{N+l}}{M}\frac{n}{2^{N+l}}-q\\ &=\frac{kn}{M2^{N+l}}+\frac nM-\frac{n-r}{M}\\ &=\frac{kn+r}{M2^{N+l}} \end{align*} \]易知 \(\dfrac{kn+r}{M2^{N+l}}>0\)
又有 \(0\le k\le 2^l,0\le n < 2^N\)，故 \(kn< 2^{N+l}\)，所以 \(\dfrac{kn+r}{M2^{N+l}}<1\)
至此已有 \(0< \dfrac{nm}{2^{N+l}}-q<1\)，所以 \(q=\lfloor\dfrac{nm}{2^{N+l}}\rfloor\)。
我们只需要一开始置 \(l\leftarrow\lfloor\log_2M\rfloor\)，之后不断缩小 \(l\)，使得只有唯一的 \(m\) 满足性质即可。
现在还有一个问题是 \(m\in[0,2^{N+l})\)，可能会超出我们的值域 \([0,2^N)\)，不过我们有
\[2^{N+l}+k=mM,k\in[0,2^l] \]显然 \(M\ge 2\)，于是就有 \(m\le\dfrac{2^{N+l}+2^l}{M}< 2^{N+1}\)，于是可以把 \(nm\) 变成 \(n(m-2^N)+n2^N\)。
至此，我们完成了整除，取模调用一下整除就好了。
还有一个细节是如果常数是偶数，我们要把所有 \(2\) 都提取出来做，这样保证了上面这个过程中的 \(M\) 是奇数。
给出奇怪的实现（怎么比不优化的慢啊）
class fastmod
{
private:
ull m;
u32 l;
inline void init()
{
if(!p) cerr<<("p is zero!\n");
l=32u-__builtin_clz(p-1);
ull low=(ull(1)<<(32u+l))/p;
ull high=((1ull<<(32u+l))+(1ull<<l))/p;
while( (low>>1)<(high>>1) && l>0 ) low>>=1,high>>=1,--l;
m=high;
}
inline u32 mulhigh(u32 x,u32 y) { return ((ull)x*y)>>32ull; }
public:
u32 p;
inline u32 div(u32 x)
{
if(m>UINT32_MAX)
{
u32 t=mulhigh(x,m-(1ull<<32ull));
return (((x-t)>>1)+t)>>(l-1);
} return mulhigh(x,m)>>l;
} inline u32 mod(u32 x) { return x-p*div(x); }
fastmod():p(998244353) { init(); }
fastmod(u32 x):p(x) { init(); }
};
class safeFastmod
{
private:
fastmod p;
u32 err;
bool al2;
inline void init()
{
if(!err) cerr<<("p is zero!\n");
al2=0;
u32 r=0;
while(!(err&1u)) err>>=1u,++r;
if(err==1u)
{
err=r,al2=1;
return;
}
p=fastmod(err);
err=r;
}
public:
inline u32 div(u32 x)
{
if(al2) return x>>err;
return p.div(x>>err);
}
inline u32 mod(u32 x)
{
if(al2) return x&((1u<<err)-1u);
return x-((div(x)*p.p)<<err);
}
inline u32 mod_num() { return al2?1u<<err:p.p<<err; }
safeFastmod() { p=fastmod(998244353),err=1,al2=0; }
safeFastmod(u32 p) { err=p; init(); }
};