• 快速幂 蒙格马利算法


    快速求幂 

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

    未完待续,,,

     

     

     

  • 相关阅读:
    通过5G网络释放触觉互联网的力量
    架构师的主要职责和一些误区
    Codeforces 305E Playing with String
    hdu3949:XOR
    bzoj1923: [Sdoi2010]外星千足虫
    bzoj1013: [JSOI2008]球形空间产生器sphere
    bzoj1770: [Usaco2009 Nov]lights 灯
    一些还没有写的AC自动机题
    bzoj2553: [BeiJing2011]禁忌
    bzoj1030: [JSOI2007]文本生成器
  • 原文地址:https://www.cnblogs.com/jeff-wgc/p/4486165.html
Copyright © 2020-2023  润新知