Problem

快速求解大数模乘。

$$ a\cdot b \mod C$$

$a, b, C$ 都很大，十进制 $n$ 位数字。其中 $a,b< C$，$C$ 为奇数。

Solution

Solution 1

我这里抛砖引玉地随便给几个我能想到的解法。

一个是先算出 $a*b$ ，然后再算出 $a*b \mod C$ 。

复杂度应该至少是 $\Theta (n^2)$。

Solution 2

快速乘。

将 b 分拆成 2 进制数组，以类似快速幂的方式进行处理。然后高精度取膜大概就是高精度整除法和高精度乘法。这俩一起优化（FFT/Karatsuba）（我不会）的话就能到 $\Theta (n^{1.5})$ 的亚子。

Montgomery Algorithm

低效的部分在哪里呢？取膜和除法。将这两个部分几乎能完全删去的算法就是蒙哥马利算法。更高效的模乘算法。

考虑将 $b$ 拆成 $k$ 位 $r$ 进制数（这里不妨假设 $r=2$ ，事实上 $r$ 一般为 $2^v$，因为这样处理更快）：
$$b = \sum_{i=0}^{k-1} b_i\cdot 2^i$$

那么则有：
$$a\cdot b = \sum_{i=0}^{k-1} a\cdot b_i\cdot 2^i = (\cdots((a\cdot b_{k-1})\cdot 2 + a\cdot b_{k-2})\cdot 2 + \cdots )\cdot 2 + a\cdot b_0$$

一个类似秦久韶算法的式子。$a\cdot b_{k-1}$ 为高精乘低精。再乘个 $r$ 即移位。所以最难的部分还是在取膜。

蒙哥马利算法提供的解在于，其不计算 $a\cdot b \mod C$ ，而选择计算 $a\cdot b \cdot r^{-k} \mod C$。

所以原来那个秦久韶的式子就变成了：
$$(\cdots((a\cdot b_{0})\cdot 2^{-1} + a\cdot b_{1})\cdot 2^{-1} + \cdots )\cdot 2^{-1}$$

所以事实上，我们每次需要计算的式子变成了：
$$(d+a\cdot b_{i})\cdot 2^{-1} \mod C$$

$d$ 是上一次迭代的值。为了简化 $X\cdot 2^{-1}$ 的操作，我们可以尝试将 $d+a\cdot b_{i}$ 加上一个 $q\cdot C$，使得 $d+a\cdot b_{i} + q\cdot C$ 能被 $2$ 整除，同时不改变原式子的值，这样乘上 $2^{-1}$ 就变成了移位的操作。现在的问题变成了如何计算 $q\cdot C$ 这个值。

令 $z=d+a\cdot b_{i}$，则我们要找到一个值 $q$ 使得 $z+q\cdot C \equiv 0 \mod r$ 。

解得 $q = -z\cdot C^{-1} \mod r$。

即 $q = (d_0+a_0\cdot b_i)\cdot -C^{-1} \mod r$。

因此，我们需要计算 $C$ 在膜 $r$ 意义下的逆元 $C^{-1}$。根据欧拉定理，$C^{-1}\equiv C^{\phi (r)-1} \mod r$，又 $\phi (r) = \phi (2^v) = 2^{v-1}$，故此处我们只需要算一下低精快速幂和低精乘法就行哒。

至此，$(d+a\cdot b_{i})\cdot 2^{-1}$ 这个式子我们已经基本搞定。$\mod C$ 就可以直接去掉了。

下证 $d<2C$ 。

假设 $d<2C$ 。而 $a< C$，且 $b_i, q\leq r-1$，故 $d+a\cdot b_{i}+q\cdot C<2C+2(r-1)C=2rC$，故新 $d'=d\cdot r^{-1}<2C$。初始时 $d=0<2C$，故根据数学归纳法，结论成立。

令函数montmul(a, b)计算 $a\cdot b \cdot r^{-k} \mod C$，令 $p = r^{2k} \mod C$。概括为以下四个步骤：

a' = montmul(a, p); // a*r^k 进入蒙哥马利域
b' = montmul(b, p); // b*r^k
res = montmul(a, b); // a*b*r^k
ans = montmul(res, 1); // a*b 退出蒙哥马利域

$p = r^{2k} = 2^{2vk}\mod C$ 使用乘法与减法暴力求幂即可。

算法复杂度约为 $\Theta (log_2^2C)$。从这个复杂度可以看出事实上，蒙哥马利算法似乎并不比 Solution 1 更优，但实际上算法耗时的大头在于计算 $p$ ，算法主要是将除法变成了移位操作，所以全部体现在常数优化得到的加速上....

蒙哥马利算法是 RSA 算法的核心。为什么要提到 RSA 呢？因为最近在做些尝试。因为除了 RSA 需要这样的底层常数优化以外也没必要提到蒙哥马利算法...

DLC - 蒙哥马利快速模幂

就是把快速幂里的暴力乘和暴力膜变成蒙哥马利快速模乘啦。

Code

此处演示无高精度，$r=2$ 的情况。

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;

int count_bit(ll x) {
    int ret=0;
    while(x) ret++, x>>=1;
    return ret;
}

int k;
ll montmul(ll a, ll b, ll n) {
    ll d=0, q;
    for(int i=0; i<k; i++) {
        q=((a&1)*((b>>i)&1)+(d&1))&1;
        d=(a*((b>>i)&1)+d+q*n)>>1;
    }
    if(d>=n) d-=n;
    return d;
}


int main() {
    ios::sync_with_stdio(false), cin.tie(0);
    ll a, x, mod, p;
    cin>>a>>x>>mod;
    ll _a=a;
    ll _x=x;
    k = count_bit(mod)+1;
    p = (1ll<<(2*k))%mod;
    ll ans = montmul(1, p, mod);
    a = montmul(a, p, mod);
    while(x) {
        if(x&1) ans=montmul(ans, a, mod);
        a=montmul(a, a, mod);
        x>>=1;
    }
    cout<<_a<<'^'<<_x<<" mod "<<mod<<'='<<montmul(ans, 1, mod);
    return 0;
}

Something else

好久没搞算法了...

如果上文有任何错误，我丝毫不会意外...

Reference / 扩展阅读

如何高效进行模乘、模幂运算？——蒙哥马利算法（Montgomery Algorithm）从入门到精通
 牛客国庆 Day5B(蒙哥马利算法) - qkoqhh
蒙哥马利算法详解 - 轩辕 233

本文链接：https://pst.iorinn.moe/archives/montgomery.html
许可： https://pst.iorinn.moe/license.html
若无特别说明，博客内的文章默认将采用 CC BY 4.0 许可协议进行许可☆

网站

加入讨论吧...