快速幂原理(二进制及二分思想)
把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,都能用二进制来表示。比如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;
}
里面的乘号,是矩阵乘的运算,res是结果矩阵。
第3行代码每进行一次,二进制数就少了最后面的一个1。二进制数有多少个1就第3行代码就执行多少次。
我们在做很多”连续“的问题的时候都会用到二进制将他们离散简化
1.多重背包问题
2.树状数组
3.状态压缩DP
……………还有很多。。。究其根本还是那句话:化连续为离散。。很多时候我们并不是为了解决一个问题而使用二进制,更多是时候是为了优化而使用它。所以如果你想让你的程序更加能适应大数据的情况,那么学习学习二进制及其算法思想将会对你有很大帮助。
详见:点击打开链接
快速幂取膜
long long qulpow(long long x,long long y,long long MOD)
{
long long tmp=1;
x=x % MOD;
while(y!=0)
{
if(y&1)
tmp=(tmp*x) % MOD;
x=(x*x) % MOD;
y>>=1;
}
return tmp;
}
矩阵相乘
const int MOD=10000;
struct mat
{
int a[2][2];//这里数据范围就用小的示范
};
mat mat_mul(mat x,mat y)//实现两个矩阵相乘,返回的还是一个矩阵。
{
mat res;//用来表示得到的新的矩阵;
memset(res.a,0,sizeof(res.a));
for(int i=0;i<2;i++)
for(int j=0;j<2;j++)
for(int k=0;k<2;k++)
{
res.a[i][j]+=x.a[i][k]*y.a[k][j];
res.a[i][j]%=MOD;//这一步看题目具体需要了
}
return res;
}
或
const int N=100;
int c[N][N];
void multi(int a[][N],int b[][N],int n)
{
memset(c,0,sizeof c);
for(int i=1;i<=n;i++)
for(int j=1;j<=n;j++)
for(int k=1;k<=n;k++)
c[i][j]+=a[i][k]*b[k][j];
}
矩阵快速幂
const long long MOD=1e9+7;
long long n,res;
struct Node
{
long long c[2][2];
}t;
Node mult(Node a, Node b)
{
Node c = {0};
for(int i = 0; i < 2; i++)
for(int j = 0; j < 2; j++)
for(int k = 0; k < 2; k++)
{
c.c[i][j] += (a.c[i][k] * b.c[k][j]) % MOD;//MOD根据题目需要选取
c.c[i][j] %= MOD;
}
return c;
}
Node pow(long long n)
{
Node pt = t;
if(n < 0) return pt;
while(n)
{
if( n & 1 )
{
pt = mult(pt, t);
n--;
}
t = mult(t, t);
n >>= 1;
}
return pt;
}
long long solve(long long n)
{
t.c[0][0] = 1;
t.c[0][1] = 1;
t.c[1][0] = 1;
t.c[1][1] = 0;
Node ans = pow(n-2);
res=ans.c[0][0] * 1;
return res;//res即为第n项斐波拉契数的结果
}
或
const int N=10;
int tmp[N][N];
void multi(int a[][N],int b[][N],int n)
{
memset(tmp,0,sizeof tmp);
for(int i=1;i<=n;i++)
for(int j=1;j<=n;j++)
for(int k=1;k<=n;k++)
tmp[i][j]+=a[i][k]*b[k][j];
for(int i=1;i<=n;i++)
for(int j=1;j<=n;j++)
a[i][j]=tmp[i][j];
}
int res[N][N];
void Pow(int a[][N],int n)
{
memset(res,0,sizeof res);
for(int i=1;i<=n;i++) res[i][i]=1;
while(n)
{
if(n&1)
multi(res,a,n);//res=res*a;复制直接在multi里面实现了;
multi(a,a,n);//a=a*a
n>>=1;
}
}
先通过入门的题目来讲应用矩阵快速幂的套路(会这题的也可以看一下套路):
例一:http://poj.org/problem?id=3070
题目:斐波那契数列f(n),给一个n,求f(n)%10000,n<=1e9;
(这题是可以用fibo的循环节去做的,不过这里不讲,反正是水题)
矩阵快速幂是用来求解递推式的,所以第一步先要列出递推式:
f(n)=f(n-1)+f(n-2)
第二步是建立矩阵递推式,找到转移矩阵:
,简写成T * A(n-1)=A(n),T矩阵就是那个2*2的常数矩阵,而
这里就是个矩阵乘法等式左边:1*f(n-1)+1*f(n-2)=f(n);1*f(n-1)+0*f(n-2)=f(n-1);
这里还是说一下构建矩阵递推的大致套路,一般An与A(n-1)都是按照原始递推式来构建的,当然可以先猜一个An,主要是利用矩阵乘法凑出矩阵T,第一行一般就是递推式,后面的行就是不需要的项就让与其的相乘系数为0。矩阵T就叫做转移矩阵(一定要是常数矩阵),它能把A(n-1)转移到A(n);然后这就是个等比数列,直接写出通项:此处A1叫初始矩阵。所以用一下矩阵快速幂然后乘上初始矩阵就能得到An,这里An就两个元素(两个位置),根据自己设置的A(n)对应位置就是对应的值,按照上面矩阵快速幂写法,res[1][1]=f(n)就是我们要求的。
给一些简单的递推式
1.f(n)=a*f(n-1)+b*f(n-2)+c;(a,b,c是常数)
2.f(n)=c^n-f(n-1) ;(c是常数)
继续例题二:poj3233
例三:HDU2276
例四:HDU 5015
当然矩阵还有很多奇葩的递推,比如这矩阵优化的dp题。
推荐一些题目:
简单的:
http://acm.hdu.edu.cn/showproblem.php?pid=1757
不简单的:
http://acm.hdu.edu.cn/showproblem.php?pid=3483
http://acm.hdu.edu.cn/showproblem.php?pid=2855
http://acm.hdu.edu.cn/showproblem.php?pid=3658
http://acm.hdu.edu.cn/showproblem.php?pid=4565