快速幂 蒙格马利算法

快速求幂 

约定:
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;
}

未完待续,,,

 

 

 

原文地址:https://www.cnblogs.com/jeff-wgc/p/4486165.html