快速幂 蒙格马利算法
快速求幂
约定：
x%y为x取模y，即x除以y所得的余数，当x<y时，x%y=x，所有取模的运算对象都为整数。
x^y表示x的y次方。
乘方运算的优先级高于乘除和取模，加减的优先级最低。
见到x^y/z这样，就先算乘方，再算除法。
A/B，称为A除以B，也称为B除A。
若A%B=0，即称为A可以被B整除，也称B可以整除A。
A*B表示A乘以B或称A乘B，B乘A，B乘以A……都一样，复习一下小学数学
公因数：两个不同的自然数A和B，若有自然数C可以整除A也可以整除B，那么C就是A和B的公因数。
公倍数：两个不同的自然数A和B，若有自然数C可以被A整除也可以被B整除，那么C就是A和B的公倍数。
互质数：两个不同的自然数，它们只有一个公因数1，则称它们互质。
费马小定理：
有N为任意正整数，P为素数，且N不能被P整除（显然N和P互质），则有：
N^P%P=N(即：N的P次方除以P的余数是N)
但是我查了很多资料见到的公式都是这个样子：
(N^(P-1))%P=1
后来分析了一下，两个式子其实是一样的，可以互相变形得到，原式可化为：
(N^P-N)%P=0(即：N的P次方减N可以被P整除，因为由费马小定理知道N的P次方除以P的余数是N)
把N提出来一个，N^P就成了你N*(N^(P-1))，那么(N^P-N)%P=0可化为：(N*(N^(P-1)-1))%P=0
请注意上式，含义是：N*(N^(P-1)-1)可以被P整除
又因为N*(N^(P-1)-1)必能整除N（这不费话么!）
所以，N*(N^(P-1)-1)是N和P的公倍数，小学知识了^_^
又因为前提是N与P互质，而互质数的最小公倍数为它们的乘积，所以一定存在正整数M使得等式成立：
N*(N^(P-1)-1)=M*N*P
两边约去N，化简之：
N^(P-1)-1=M*P
因为M是整数，显然：
(N^(P-1)-1)%P=0
即：
N^(P-1)%P=1
============================================
积模分解公式
先有一个引理，如果有：X%Z=0，即X能被Z整除，则有：
(X+Y)%Z=Y%Z
这个不用证了吧...
设有X、Y和Z三个正整数，则必有：(X*Y)%Z=((X%Z)*(Y%Z))%Z
想了很长时间才证出来，要分情况讨论才行：
1.当X和Y都比Z大时，必有整数A和B使下面的等式成立：
X=Z*I+A（1）
Y=Z*J+B（2）
不用多说了吧，这是除模运算的性质！
将（1）和（2）代入(X*Y)modZ得：((Z*I+A)(Z*J+B))%Z
乘开，再把前三项的Z提一个出来，变形为：(Z*(Z*I*J+I*A+I*B)+A*B)%Z（3）
因为Z*(Z*I*J+I*A+I*B)是Z的整数倍……晕，又来了。
概据引理，（3）式可化简为：(A*B)%Z
又因为：A=X%Z，B=Y%Z，代入上面的式子，就成了原式了。
2.当X比Z大而Y比Z小时，一样的转化：
X=Z*I+A
代入(X*Y)%Z得：
(Z*I*Y+A*Y)%Z
根据引理，转化得：(A*Y)%Z
因为A=X%Z，又因为Y=Y%Z，代入上式，即得到原式。
同理，当X比Z小而Y比Z大时，原式也成立。
3.当X比Z小，且Y也比Z小时，X=X%Z，Y=Y%Z，所以原式成立。
=====================================================
快速计算乘方的算法
如计算2^13，则传统做法需要进行12次乘法。
/*计算n^p*/
unsigned power(unsigned n,unsigned p)
{
for(int i=0;i<p;i++) n*=n;
return n;
}
该死的乘法，是时候优化一下了！把2*2的结果保存起来看看，是不是成了：4*4*4*4*4*4*2
再把4*4的结果保存起来：16*16*16*2
一共5次运算，分别是2*2、4*4和16*16*16*2
这样分析，我们算法因该是只需要计算一半都不到的乘法了。
为了讲清这个算法，再举一个例子2^7：2*2*2*2*2*2*2
两两分开：(2*2)*(2*2)*(2*2)*2
如果用2*2来计算，那么指数就可以除以2了，不过剩了一个，稍后再单独乘上它。
再次两两分开，指数除以2： ((2*2)*(2*2))*(2*2)*2
实际上最后一个括号里的2 * 2是这回又剩下的，那么，稍后再单独乘上它
现在指数已经为1了，可以计算最终结果了：16*4*2=128
优化后的算法如下：
unsigned Power(unsigned n,unsigned p)
{
unsigned main=n; //用main保存结果
unsigned odd=1; //odd用来计算“剩下的”乘积
while (p>1)
{//一直计算，直到指数小于或等于1
if((p%2)!=0)
{// 如果指数p是奇数，则说明计算后会剩一个多余的数，那么在这里把它乘到结果中
odd*=main; //把“剩下的”乘起来
}
main*=main; //主体乘方
p/=2; //指数除以2
}
return main*odd; //最后把主体和“剩下的”乘起来作为结果
}
够完美了吗？不，还不够！看出来了吗？main是没有必要的，并且我们可以有更快的代码来判断奇数。要知道除法或取模运算的效率很低，所以我们可以利用偶数的一个性质来优化代码，那就是偶数的二进制表示法中的最低位一定为0!
完美版：
unsigned Power(unsigned n, unsigned p)
{ // 计算n的p次方
unsigned odd = 1; //oddk用来计算“剩下的”乘积
while (p > 1)
{ // 一直计算到指数小于或等于1
if (( p & 1 )!=0)
{ // 判断p是否奇数，偶数的最低位必为0
odd *= n; // 若odd为奇数，则把“剩下的”乘起来
}
n *= n; // 主体乘方
p /= 2; // 指数除以2
}
return n * odd; // 最后把主体和“剩下的”乘起来作为结果
}
========================================================
蒙格马利”快速幂模算法
后面我们会用到这样一种运算：(X^Y)%Z
问题是当X和Y很大时，只有32位的整型变量如何能够有效的计算出结果？
考虑上面那份最终的优化代码和再上面提到过的积模分解公式，我想你也许会猛拍一下脑门，吸口气说：“哦，我懂了！”。
下面的讲解是给尚没有做出这样动作的同学们准备的。X^Y可以看作Y个X相乘，即然有积模分解公式，那么我们就可以把Y个X相乘再取模的过程分解开来，比如：(17^25)%29则可分解为：( ( 17 * 17 ) % 29 * ( 17 * 17 ) % 29 * ……
如果用上面的代码将这个过程优化，那么我们就得到了著名的“蒙格马利”快速幂模算法：
unsigned Montgomery(unsigned n, unsigned p, unsigned m)
{ // 快速计算 (n ^ e) % m 的值，与power算法极类似
unsigned r = n % m; // 这里的r可不能省
unsigned k = 1;
while (p > 1)
{
if ((p & 1)!=0)
{
k = (k * r) % m; // 直接取模
}
r = (r * r) % m; // 同上
p /= 2;
}
return (r * k) % m; // 还是同上
}
上面的代码还可以优化。下面是蒙格马利极速版：
unsigned Montgomery(unsigned n,unsigned p,unsigned m)
{ //快速计算(n^e)%m的值
unsignedk=1;
n%=m;
while(p!=1)
{
if(0!=(p&1))k=(k*n)%m;
n=(n*n)%m;
p>>=1;
}
return(n*k)%m;
}
未完待续，，，
