模重复平方算法

算法原理

模重复平方算法是用来快速计算 b n m o d    m b^n \mod m bnmodm的一个算法。
考虑直接计算 b n m o d    m b^n \mod m bnmodm，需要 n − 1 n-1 n−1次乘法，也就是递归计算 b n ≡ ( b n − 1 m o d    m ) ⋅ b m o d    m . b^n \equiv (b^{n-1} \mod m) \cdot b \mod m. bn≡(bn−1modm)⋅bmodm.不过，当 n n n很大的时候，计算会非常耗时。
现在，考虑将 n n n的二进制表示 n = n 0 + n 1 2 + n 2 2 2 + ⋯ + n k − 1 2 k − 1 n=n_0+n_12+n_22^2+\cdots +n_{k-1}2^{k-1} n=n0​+n1​2+n2​22+⋯+nk−1​2k−1，其中 n 0 , n 1 , ⋯   , n k − 1 n_0,n_1,\cdots,n_{k-1} n0​,n1​,⋯,nk−1​为0或者1.
则:
b n ≡ b n 0 + n 1 2 + n 2 2 2 + ⋯ + n k − 1 2 k − 1 m o d    m ≡ b n 0 ⋅ ( b 2 ) n 1 ⋅ ( ( b 2 ) 2 ) n 2 ⋯ ( b 2 k − 1 ) n k − 1 m o d    m \begin{aligned} b^n& \equiv b^{n_0+n_12+n_22^2+\cdots +n_{k-1}2^{k-1}} \mod m\ &\equiv b^{n_0} \cdot (b^2)^{n_1}\cdot ((b^2)^2)^{n_2} \cdots (b^{2^{k-1}})^{n_{k-1}} \mod m \end{aligned} bn​≡bn0​+n1​2+n2​22+⋯+nk−1​2k−1modm≡bn0​⋅(b2)n1​⋅((b2)2)n2​⋯(b2k−1)nk−1​modm​
通过观察，我们可以知道，我们需要 b b b的偶次幂。设 b 0 = b b_0=b b0​=b的话，下一个偶次幂 b 1 = b 0 2 m o d    m b_1=b_0^2 \mod m b1​=b02​modm，以此类推， b 2 = b 1 2 m o d    m , ⋯   , b k − 1 = b k − 2 2 m o d    m b_2=b_1^2 \mod m, \cdots ,b_{k-1}=b_{k-2}^2 \mod m b2​=b12​modm,⋯,bk−1​=bk−22​modm，然后代入上面的公式，有：
b n ≡ b 0 n 0 ⋅ b 1 n 1 ⋅ b 2 n 2 ⋯ b k − 1 n k − 1 m o d    m . b^n \equiv b_0^{n_0} \cdot b_1^{n_1} \cdot b_2^{n_2} \cdots b_{k-1}^{n_{k-1}} \mod m. bn≡b0n0​​⋅b1n1​​⋅b2n2​​⋯bk−1nk−1​​modm.
又由于 n 0 , n 1 , ⋯   , n k − 1 n_0,n_1,\cdots,n_{k-1} n0​,n1​,⋯,nk−1​为0或者1,所以，我们只需要通过不断平方计算出 b i b_i bi​，令初始结果为1，然后，如果 n i n_i ni​为1则乘上 b i m o d    m b_i \mod m bi​modm，否则，不管，继续计算 b i + 1 b_{i+1} bi+1​.
可以看出，通过上面的计算，我们只需要最多 2 [ l o g 2 n ] 2[log_2n] 2[log2​n]次乘法，和一次 n n n的分解。

c++实现代码

#include
using namespace std;
int msa(int b,int n,int m){
	int b1=b;
	int res=1;
	while(n>0){
		if(n%2==1){
			res=res*b1%m;
		}
		b1=b1*b1%m;
		n=n/2;
	}
	return res;
}
int main(){
	int b,n,m;
	cout<<"Input b:"<<endl;
	cin>>b;
	cout<<"Input n:"<<endl;
	cin>>n;
	cout<<"Input m:"<<endl;
	cin>>m;
	cout<<"The result of b^n mod m is:"<<endl;
	cout<<msa(b,n,m)<<endl;
}

gmp函数

在gmp库中，实现了大整数的幂运算，其函数原型如下：
void mpz_powm (mpz_t rop, const mpz_t base, const mpz_t exp, const mpz_t mod)
rop是返回的结果。