算法原理

模重复平方算法是用来快速计算$b^n\mathrm{m}\mathrm{o}\mathrm{d}m$的一个算法。
考虑直接计算$b^n\mathrm{m}\mathrm{o}\mathrm{d}m$，需要$n-1$次乘法，也就是递归计算$$b^n\equiv (b^{n-1}\mathrm{m}\mathrm{o}\mathrm{d}m)\cdot b\mathrm{m}\mathrm{o}\mathrm{d}m.$$不过，当$n$很大的时候，计算会非常耗时。
现在，考虑将$n$的二进制表示$n=n_0+n_{1}2+n_2{2}^2+\cdots +n_{k-1}{2}^{k-1}$，其中$n_0,n_1,\cdots ,n_{k-1}$为0或者1.
则:
$\begin{array}{rl}{b}^n& \equiv b^{n_0+n_{1}2+n_2{2}^2+\cdots +n_{k-1}{2}^{k-1}}\mathrm{m}\mathrm{o}\mathrm{d}m\\ & \equiv b^{n_0}\cdot (b^2{)}^{n_1}\cdot ((b^2{)}^2{)}^{n_2}\cdots (b^{2^{k-1}}{)}^{n_{k-1}}\mathrm{m}\mathrm{o}\mathrm{d}m\end{array}$
通过观察，我们可以知道，我们需要$b$的偶次幂。设$b_0=b$的话，下一个偶次幂$b_1=b_0^2\mathrm{m}\mathrm{o}\mathrm{d}m$，以此类推，$b_2=b_1^2\mathrm{m}\mathrm{o}\mathrm{d}m,\cdots ,b_{k-1}=b_{k-2}^2\mathrm{m}\mathrm{o}\mathrm{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}}\mathrm{m}\mathrm{o}\mathrm{d}m.$$
又由于$n_0,n_1,\cdots ,n_{k-1}$为0或者1,所以，我们只需要通过不断平方计算出$b_i$，令初始结果为1，然后，如果$n_i$为1则乘上$b_i\mathrm{m}\mathrm{o}\mathrm{d}m$，否则，不管，继续计算$b_{i+1}$.
可以看出，通过上面的计算，我们只需要最多$2[lo{g}_{2}n]$次乘法，和一次$n$的分解。

c++实现代码

#include<iostream>
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是返回的结果。