取模从古以来就是一个难题。
比如有一个数 \(p\),和两个数 \(0\leq x,y<p\),计算 \((xy\bmod p)\)。暴力做除法的复杂度是 \(O(\dfrac{\log x\log p}{w})\),这里请注意。高精度乘除的复杂度是不同的,这是因为乘法可以把两边都压位,而除法必须枚举一个。接下来,我们介绍 Montgomery 形式,不仅可以优化取模,还能干其它事。
有函数 \(R(x)\),定义域是 \([0,p)\),返回 \(xr\bmod p\)。\(R'(x)\),定义域是 \([0,p^2)\),返回 \(\dfrac{x}{r}\bmod p\),注意它们的值域都是 \([0,p)\),可是定义域不同,因此它们并不互为反函数。
如何求 \(R'(x)\)?
答:设 \(\dfrac{x}{r}\equiv t\bmod p\),则 \(x\equiv tr\bmod p\)。改写为 \(x=tr+kp\)。你现在要求出一个 \([0,p)\) 的 \(t\)。
换模数(这也是数论问题常见的处理技巧?),改为两边同时 \(\bmod r\),则 \(kp\equiv x\bmod r\),解得 \(k=(xp^{-1}\bmod r)\)。
这样我们得到了一个 \([0,r)\) 的 \(k\)。又知 \(t=\dfrac{x-kp}{r}\),\(x\in [0,p^2)\),所以有 \(|t|<\dfrac{p^2+rp}{r}\)。
请注意,刚才我们并没有说 \(r\) 是几。但是从中可以看出特点:要好求逆(也说明要和 \(p\) 互素,这里 \(p\) 必须是质数),且 \(r>p\)。最好的办法是让 \(r\) 取成一个 \(2^{2^q}\) 形式。则 \(-p<t<2p\),稍作修改即可将 \(t\) reduce 为 \([0,p)\) 的数。取模关键部分在这里。
下一个问题:如何求 \(R(x)\)?
因为 \(x\equiv \dfrac{xr^2}{r}\bmod p\),所以改为算 \(R'(xz)\) 即可,其中 \(z=(r^2\bmod p)\),可以预处理!
请注意,下面要说一个关键点。我们从此考虑那些运算时,不要 再把那些原来的 \(x\) 带入计算了,通通改成 \(R(x)\)。这为什么是可行的?
第一点:\(R(x\pm y)=R(x)\pm R(y)\),reduce 是极其容易的。
第二点:\(R(xy)\equiv\dfrac{R(x)R(y)}{r}\bmod p\),因此有 \(R(xy)=R'(R(x)R(y))\)。这里注意 \(R'\) 的定义域可以有平方,因此直接一个乘法一个 reduce 就好了。
第三点:也就是最后怎么还原?其实也非常容易:\(x=R'(R(x))\)。
于是就赢麻了。假如我们高精度都是在 \(2\) 进制意义下计算,则对一个 \(2\) 的次幂取模是轻松的(大概相当于 resize)。最后一个问题,计算一个奇数在 \(\bmod 2^{2^q}\) 意义下的逆元。
这东西像不像多项式求逆?没错,牛顿迭代即可。
设要求出 \(a\) 的逆元,已经有的数是 \(b\),每轮 \(b\leftarrow 2b-b^2a\) 即可。因为倍增,所以这部分复杂度 \(O(\dfrac{\log^2 p}{w^2})\),很可能不足以成为复杂度瓶颈。
每次取模的复杂度降为了 \(O(\dfrac{\log^2 p}{w^2}\))。