快速幂取模算法

在网站上一直没有找到有关于快速幂算法的一个详细的描述和解释，这里，我给出快速幂算法的完整解释，用的是C语言，不同语言的读者只好换个位啦，毕竟读C的人较多~

所谓的快速幂，实际上是快速幂取模的缩写，简单的说，就是快速的求一个幂式的模(余)。在程序设计过程中，经常要去求一些大数对于某个数的余数，为了得到更快、计算范围更大的算法，产生了快速幂取模算法。[有读者反映在讲快速幂部分时有点含糊，所以在这里对本文进行了修改，作了更详细的补充，争取让更多的读者一目了然]

我们先从简单的例子入手：求= 几。

算法1.首先直接地来设计这个算法：

int ans = 1;

for(int i = 1;i<=b;i++)

{

ans = ans * a;

}

ans = ans % c;

这个算法的时间复杂度体现在for循环中，为O（b）.这个算法存在着明显的问题，如果a和b过大，很容易就会溢出。

那么，我们先来看看第一个改进方案：在讲这个方案之前，要先有这样一个公式：

.这个公式大家在离散数学或者数论当中应该学过，不过这里为了方便大家的阅读，还是给出证明：

引理1：

上面公式为下面公式的引理，即积的取余等于取余的积的取余。

证明了以上的公式以后，我们可以先让a关于c取余，这样可以大大减少a的大小，

于是不用思考的进行了改进：

算法2：

int ans = 1;

a = a % c; //加上这一句

for(int i = 1;i<=b;i++)

{

ans = ans * a;

}

ans = ans % c;

聪明的读者应该可以想到，既然某个因子取余之后相乘再取余保持余数不变，那么新算得的ans也可以进行取余，所以得到比较良好的改进版本。

算法3：

int ans = 1;

a = a % c; //加上这一句

for(int i = 1;i<=b;i++)

{

ans = (ans * a) % c;//这里再取了一次余

 

}

ans = ans % c;

这个算法在时间复杂度上没有改进，仍为O(b)，不过已经好很多的，但是在c过大的条件下，还是很有可能超时，所以，我们推出以下的快速幂算法。

快速幂算法依赖于以下明显的公式，我就不证明了。

有了上述两个公式后，我们可以得出以下的结论：

1.如果b是偶数，我们可以记k = a2 mod c，那么求(k)b/2 mod c就可以了。

2.如果b是奇数，我们也可以记k = a2 mod c，那么求

((k)b/2 mod c × a ) mo d c =((k)b/2 mod c * a) mod c 就可以了。

 

那么我们可以得到以下算法：

算法4：

int ans = 1;

a = a % c;

if(b%2==1)

ans = (ans * a) mod c; //如果是奇数，要多求一步，可以提前算到ans中

k = (a*a) % c; //我们取a2而不是a

for(int i = 1;i<=b/2;i++)

{

ans = (ans * k) % c;

}

ans = ans % c;

 

我们可以看到，我们把时间复杂度变成了O(b/2).当然，这样子治标不治本。但我们可以看到，当我们令k = (a * a) mod c时，状态已经发生了变化，我们所要求的最终结果即为(k)b/2 mod c而不是原来的ab mod c，所以我们发现这个过程是可以迭代下去的。当然，对于奇数的情形会多出一项a mod c，所以为了完成迭代，当b是奇数时，我们通过

ans = (ans * a) % c;来弥补多出来的这一项，此时剩余的部分就可以进行迭代了。

 

形如上式的迭代下去后，当b=0时，所有的因子都已经相乘，算法结束。于是便可以在O（log b）的时间内完成了。于是，有了最终的算法：快速幂算法。

算法5：快速幂算法

 

int ans = 1;

a = a % c;

while(b>0)

{

 

if(b % 2 == 1)

ans = (ans * a) % c;

b = b/2;

a = (a * a) % c;

}

将上述的代码结构化，也就是写成函数：

int PowerMod(int a, int b, int c)

{

int ans = 1;

a = a % c;

while(b>0)

{

 

if(b % 2 = = 1)

ans = (ans * a) % c;

b = b/2;

a = (a * a) % c;

}

return ans;

}

本算法的时间复杂度为O（logb），能在几乎所有的程序设计（竞赛）过程中通过，是目前最常用的算法之一。

以下内容仅供参考：

扩展：有关于快速幂的算法的推导，还可以从另一个角度来想。

=？ 求解这个问题，我们也可以从进制转换来考虑：

将10进制的b转化成2进制的表达式：

那么，实际上，.

所以

注意此处的要么为0，要么为1，如果某一项，那么这一项就是1，这个对应了上面算法过程中b是偶数的情况，为1对应了b是奇数的情况[不要搞反了，读者自己好好分析，可以联系10进制转2进制的方法]，我们从依次乘到。对于每一项的计算，计算后一项的结果时用前一项的结果的平方取余。对于要求的结果而言，为时ans不用把它乘起来，[因为这一项值为1]，为1项时要乘以此项再取余。这个算法和上面的算法在本质上是一样的，读者可以自行分析，这里我说不多说了，希望本文有助于读者掌握快速幂算法的知识点，当然，要真正的掌握，不多练习是不行的。

 

 

 

矩阵 快速求幂 blog1

原文： http://www.cnblogs.com/yan-boy/archive/2012/11/29/2795294.html

矩阵的快速幂是用来高效地计算矩阵的高次方的。将朴素的o（n）的时间复杂度，降到log（n）。

这里先对原理（主要运用了矩阵乘法的结合律）做下简单形象的介绍：

一般一个矩阵的n次方，我们会通过连乘n-1次来得到它的n次幂。

但做下简单的改进就能减少连乘的次数，方法如下：

把n个矩阵进行两两分组，比如：A*A*A*A*A*A  =>  (A*A)*(A*A)*(A*A)

这样变的好处是，你只需要计算一次A*A，然后将结果(A*A)连乘自己两次就能得到A^6，即(A*A)^3=A^6。算一下发现这次一共乘了3次，少于原来的5次。

其实大家还可以取A^3作为一个基本单位。原理都一样：利用矩阵乘法的结合律，来减少重复计算的次数。

以上都是取一个具体的数来作为最小单位的长度，这样做虽然能够改进效率，但缺陷也是很明显的，取个极限的例子（可能有点不恰当，但基本能说明问题），当n无穷大的时候，你现在所取的长度其实和1没什么区别。所以就需要我们找到一种与n增长速度”相适应“的”单位长度“，那这个长度到底怎么去取呢？？？这点是我们要思考的问题。

有了以上的知识，我们现在再来看看，到底怎么迅速地求得矩阵的N次幂。

既然要减少重复计算，那么就要充分利用现有的计算结果咯！~怎么充分利用计算结果呢？？？这里考虑二分的思想。。

大家首先要认识到这一点：任何一个整数N，都能用二进制来表示。。这点大家都应该知道，但其中的内涵真的很深很深（这点笔者感触很深，在文章的最后，我将谈谈我对的感想）！！

计算机处理的是离散的信息，都是以0,1来作为信号的处理的。可想而知二进制在计算机上起着举足轻重的地位。它能将模拟信号转化成数字信号，将原来连续的实际模型，用一个离散的算法模型来解决。  好了，扯得有点多了，不过相信这写对下面的讲解还是有用的。

回头看看矩阵的快速幂问题，我们是不是也能把它离散化呢？比如A^19  =>  （A^16）*（A^2）*（A^1），显然采取这样的方式计算时因子数将是log(n)级别的(原来的因子数是n)，不仅这样，因子间也是存在某种联系的，比如A^4能通过(A^2)*(A^2)得到，A^8又能通过(A^4)*(A^4)得到，这点也充分利用了现有的结果作为有利条件。下面举个例子进行说明：

现在要求A^156,而156(10)=10011100(2) 

也就有A^156=>(A^4)*(A^8)*(A^16)*(A^128)  考虑到因子间的联系，我们从二进制10011100中的最右端开始计算到最左端。细节就说到这，下面给核心代码：

 

 

 

while(N)

 {

                if(N&1)

                       res=res*A;

                n>>=1;

                A=A*A;

 }

 

改进版（取模）a^b mod n ：

int modexp(int a, int b, int n)

{

   int t = a, ret = 1;

   while(b != 0)

  {

     if(b % 2 == 1) ret *= t % n;

     t = t * t % n;

     b /= 2;

  } 

  return ret;

}

 

 

 

 

 

里面的乘号，是矩阵乘的运算，res是结果矩阵。

第3行代码每进行一次，二进制数就少了最后面的一个1。二进制数有多少个1就第3行代码就执行多少次。

好吧，矩阵快速幂的讲解就到这里吧。在文章我最后给出我实现快速幂的具体代码（代码以3*3的矩阵为例）。

现在我就说下我对二进制的感想吧：

我们在做很多”连续“的问题的时候都会用到二进制将他们离散简化

1.多重背包问题

2.树状数组

3.状态压缩DP

……………还有很多。。。究其根本还是那句话：化连续为离散。。很多时候我们并不是为了解决一个问题而使用二进制，更多是时候是为了优化而使用它。所以如果你想让你的程序更加能适应大数据的情况，那么学习学习二进制及其算法思想将会对你有很大帮助。

 

 

 

 

快速幂或者矩阵快速幂 blog 2 

原文：http://blog.csdn.net/hikean/article/details/9749391

快速幂或者矩阵快速幂在算大指数次方时是很高效的，他的基本原理是二进制，下面的A可以是一个数也可以是一个矩阵（本文特指方阵），若是数就是快速幂算法，若是矩阵就是矩阵快速幂算法，用c++只需把矩阵设成一个类就可以，然后重载一下乘法就可以，注意为矩阵时则ANS=1，应该是ANS=E，E是单位矩阵，即主对角线是1其余的部分都是0的特殊方阵了。

 举个例子若你要算A^7你会怎么算一般你会用O(N)的算法A^7=A*A*A*A*A*A*A也许你觉得这并不慢但是若要你算A^10000000000000000呢，是不是会觉得O(N)的算法也太慢了吧这不得算死我啊，计算机也不想算了，因为有更高效的算法我们把A的指数写成二进制，这样就有了

A^7=A^111(2)  现在我们可以这么算 令ANS=1;MULTI=A  ,N=7

while(N)

{

if(N%2==1) //亦可以写成N&1 或N%2

ANS*=MULTI;

MULTI*=MULTI;

N/=2;//c++中可以写成 N>>=1;直接用位运算更快

}

写出上面的代码的执行过程就是

ANS=1;MULTI=A; 

N=7 ;N%2=1;   ANS*=MULTI; 所以ANS=A;  MULTI*=MULTI; 所以MULTI=A^2

然后 N/=2;N=3; N%2=1; ANS*=MULTI; 所以 ANS=A*A^2=A^3 ; 又MULTI*=MULTI; 所以MULTI=A^4

然后N/=2;N=1;N%2=1;ANS*=MULTI; 所以 ANS=A*A^2*A^4=A^7;又MULTI*=MULTI; 所以MULTI=A^8

然后N/=2;N=0;算法结束  是不是很巧妙呢，实际上用的乘法次数是 6次你可能觉得，那个A^7=A*A*A*A*A*A*A，不也是用了6次乘法吗有什么区别？

那是因为这个算法是log2(n)   (表示以2为底n的对数) 的复杂度，还有一个系数，大约是2 实际上计算次数就是 2*log2(n) 而普通的连乘计算的复杂度是n 乘法计算次数是n-1

这样在n很小时差别不大，但随着n的增长差距会迅速扩大，例如 n=1024时 普通方法得计算1023次乘法，但快速幂最多（因为当上面的程序执行时N的中间结果为偶数那么  ANS*=MULTI，将不会被执行，故实际的计算次数要小于 2*log2(n)）只算2*log2(n) =20次乘法，是不是很快！！！！！！！！！！

但是为什么呢？好像还有点不懂。。

实际上A^7=A^1*A^2*A^4这样每次计算乘法乘的因子都是递增的，而且还是指数递增，还有这些因子是可以递推产生的就是可以利用上次的计算每次平方就可以了，这中其实是使用的二进制的思想，因为任意一个数都可以，表示成二进制，故 A^N以定可以写成

A^(一个二进制数如101010)=A^(100000)*A^(00000)*A(1000)*A^(000)*A^(10)*A^(0)=A^(2^5)*A^(2^3)*A^(2^1)

而我们的MULTI 其实是一个数列 A^1，A^2，A^4，A^8，A^16........即A^(2^0),A^(2^1),A^(2^2)，A^(2^3),A^(2^4).......................注意到他的指数都是二进制的位权（不知道是不是这个名词，就像十进制的位权是 1 10 100 1000 10000，一样如1243=1*1000+2*100+3*10+4*1；而二进制的1011 是 1*2^(3)+0*2^(2)+1*2^(1)+1*2^(0) 这样是不是应该理解位权了呢？）实际上任何一个A^N都可以写成这个数列的某些项的乘积，因为N始终都可以表示成二进制，而把N表示成二进制后如果某项为1则说明需要乘上MULTI 否则不用乘上MULTI

于是就有了上面的代码，，，，哎怎么感觉我说的还是很不清楚呢?那就没办法

下面附上代码，另外一般要用快速幂的题都要取模，一般保证模的平方不超数的范围， 因为指数太大的数是会爆掉int 和long long 的

 

#include<iostream>

using namespace std;

#define mod 1000000007

long long quick_pow(int n,int base)

//n是指数 base是底 即计算的是base^n 当然结果是取模了的

{

    long long ans=1;//默认ans大于等于1因为不能算负指数

    long long multi=base;

    while(n)

    {

        if(n%2) ans*=multi;

        ans%=mod;//由于数太大一般要取模

        n/=2;

        multi*=multi;

        multi%=mod;

    }

return ans;

}

 

int main()

{

    int n,base;

    while(cin>>n>>base)

        cout<<quick_pow(n,base)<<endl;

    return 0;

}

可能你会问了这个算法有什么用呢？其实用的更多是使用矩阵快速幂，算递推式，注意是递推式，简单的如斐波那契数列的第一亿项的结果模上10000000后是多少你还能用递推式去，逐项递推吗？当然不能，这里就可以发挥矩阵快速幂的神威了，那斐波那契数列和矩阵快速幂能有一毛钱的关系？答案是有而且很大

斐波那契的定义是f(1)=f(2)=1; 然后f(n)=f(n-1)+f(n-2) (n>=2) 我们也可以这样定义f(1)=f(2)=1; [f(n),f(n-1)]=[f(n-1),f(n-2)][1,1,1,0],其中[1,1,1,0] 是一个2*2的矩阵 上面一行是1,1，下面一行是1,0，这样就可以化简了写成[f(n),f(n-1)]=[f(2),f(1)]*[1,1,1,0]^(n-2)

化简一下

这样就可以用矩阵快速幂，快速的推出斐波那契数列的第一亿项的值了（当然是取模的值了）是不是很神奇，类似的递推式也可以，化成这种形式，用矩阵快速幂进行计算

下面附一个矩阵快速幂的代码，当然所有矩阵都是要模的

 

# include<cstdio>

# include<cstring>

using namespace std;

#define NUM 50

int MAXN,n,mod;

struct Matrix//矩阵的类

{

int a[NUM][NUM];

void init()           //将其初始化为单位矩阵

{

memset(a,0,sizeof(a));

for(int i=0;i<MAXN;i++)

a[i][i]=1;

}

} A;

Matrix mul(Matrix a,Matrix b)  //(a*b)%mod  矩阵乘法

{

Matrix ans;

for(int i=0;i<MAXN;i++)

for(int j=0;j<MAXN;j++)

{

ans.a[i][j]=0;

for(int k=0;k<MAXN;k++)

ans.a[i][j]+=a.a[i][k]*b.a[k][j];

ans.a[i][j]%=mod;

}

return ans;

}

 

Matrix add(Matrix a,Matrix b)  //(a+b)%mod  //矩阵加法

{

int i,j,k;

Matrix ans;

for(i=0;i<MAXN;i++)

for(j=0;j<MAXN;j++)

{

ans.a[i][j]=a.a[i][j]+b.a[i][j];

ans.a[i][j]%=mod;

}

return ans;

}

 

Matrix pow(Matrix a,int n)    //(a^n)%mod  //矩阵快速幂

{

Matrix ans;

ans.init();

while(n)

{

if(n%2)//n&1

ans=mul(ans,a);

n/=2;

a=mul(a,a);

}

return ans;

}

 

Matrix sum(Matrix a,int n)  //(a+a^2+a^3....+a^n)%mod// 矩阵的幂和

{

int m;

Matrix ans,pre;

if(n==1)

return a;

m=n/2;

pre=sum(a,m);                      //[1,n/2]

ans=add(pre,mul(pre,pow(a,m)));   //ans=[1,n/2]+a^(n/2)*[1,n/2]

if(n&1)

ans=add(ans,pow(a,n));          //ans=ans+a^n

return ans;

}

 

void output(Matrix a)//输出矩阵

{

for(int i=0;i<MAXN;i++)

for(int j=0;j<MAXN;j++)

printf("%d%c",a.a[i][j],j==MAXN-1?'\n':' ');

}

int main()

{

freopen("in.txt","r",stdin);

Matrix ans;

scanf("%d%d%d",&MAXN,&n,&mod);

for(int i=0;i<MAXN;i++)

for(int j=0;j<MAXN;j++)

{

scanf("%d",&A.a[i][j]);

A.a[i][j]%=mod;

}

ans=sum(A,n);

output(ans);

return 0;

}

 

 

快速幂乘：

http://blog.csdn.net/cambridgeacm/article/details/7703809