1.引入
//整数快速幂
为了引入矩阵的快速幂,我先看了求整数的幂,即x^n
例如求x^8 就是x*x*x*x*x*x*x*x
正常的运算方式是,x的值一个个往上乘上去,乘法运算运行7次
(x x)(x x) (x x)(x x)
也可以采用这种运算方式,先进行乘法得到x^2再对x^2进行三次乘法。这种运算要明显比第一种情况要快
所以对于整数快速幂,也是结合了这种思想
(x^m)*(x^n)=x^(m+n)
x^19=(x^16)(x^2)(x^1)
//整数快速幂
int QuickPow(int x,int N)
{
int res = x;
int ans = 1;
while(N)
{
if(N&1)
{
ans = ans * res;
}
res = res*res;
N = N>>1;
}
return ans;
}
我们用 x^19来走这个函数
举例:x^19
19的二进制为:1 0 0 1 1
初始:
ans = 1 ; res = x ;
则10011最后1位是1,为奇数
ans = res *ans =x ;
res =res * res = x ^ 2 ;
然后右移一位,1001
最后一位为1,为奇数
ans = res*ans = x*(x^2) = x^3;
res = res*res = x^2*x^2 = x^4 ;
然后右移一位, 100
为偶数
res = res*res = x^4*x^4 = x^8;
然后右移一位,10
为偶数
res=res*res=x^8*x^8=x^16;
然后右移一位,1
为奇数
ans = ans*res = (x^3)*(x^16) = x^19;
res = res*res = x^32
res=x^m,m始终是与二进制位置上的权值是相互对应的。
当二进制为0时候,我们只让res*res使得幂指数乘2,对应下一个二进制权值;
当二进制为1时候,我们让ans=ans*res,乘上了x幂次。
2.矩阵快速幂
看完了整数快速幂之后,再来学习一下,矩阵快速幂,
设有一个n*n的方阵A。
给出一个数M,即计算矩阵A的M次幂 A^M
struct Matrix //结构体,矩阵类型
{
int m[maxn][maxn];
}ans,res;
//计算矩阵乘法的函数,参数是矩阵a,矩阵b和一个n,代表这两个矩阵是几阶方阵
Maxtrix Mul(Maxtrix A,Maxtrix B,int n)
{
Matrix tmp;//定义一个临时的矩阵,存放A*B的结果
for(int i=1;i<=n;i++)
for(int j=1;j<=n;j++)
tmp.m[i][j]=0;
for(int i=1;i<=n;i++)
for(int j=1;j<=n;j++)
for(int k=1;k<=n;k++)
tmp.m[i][j]+=A.m[i][k]*B.m[k][j];
return tmp;
}
//快速幂算法,求矩阵res的N次幂
void QuickPower(int N,int n)
{
//对于整数求幂来说,ans初始化为1,对于矩阵乘法来说,ans应初始化为单位阵
//对于单位矩阵来说,任何矩阵A*E=A
for(int i=1;i<=n;i++)
for(int j=1;j<=n;j++)
{
if(i==j) ans.m[i][j]=1;
else ans.m[i][j]=0;
}
while(N)
{
if(N&1)
ans=Mul(ans,res);
res=Mul(res,res);
N=N>>1;
}
}
这是简单的计算矩阵幂,在实际中的运用的话,
3.实际应用
我们以最为常见的斐波那契数列为例:
f[n]=f[n-1]+f[n-2]
f[0]=0,f[1]=1;
斐波那契数列是一个递推式求出来的
第n+1项由第n项和第n-1项递推而来
因此可以用矩阵表示:
POJ3070/NYOJ148
给出了矩阵相乘的定义,要你求出斐波那契的第n项对1e4取余。
#include <iostream>
#include <cstddef>
#include <cstring>
#include <vector>
using namespace std;
typedef long long ll;
const int mod=10000;
typedef vector<ll> vec;
typedef vector<vec> mat;
mat mul(mat &a,mat &b)//表示不会这样用,,,,
{
mat c(a.size(),vec(b[0].size()));
for(int i=0; i<2; i++)
{
for(int j=0; j<2; j++)
{
for(int k=0; k<2; k++)
{
c[i][j]+=a[i][k]*b[k][j];
c[i][j]%=mod;
}
}
}
return c;
}
mat pow(mat a,ll n)
{
mat res(a.size(),vec(a.size()));
for(int i=0; i<a.size(); i++)
res[i][i]=1;//单位矩阵;
while(n)
{
if(n&1) res=mul(res,a);
a=mul(a,a);
n/=2;
}
return res;
}
ll solve(ll n)
{
mat a(2,vec(2));
a[0][0]=1;
a[0][1]=1;
a[1][0]=1;
a[1][1]=0;
a=pow(a,n);
return a[0][1];//也可以是a[1][0];
}
int main()
{
ll n;
while(cin>>n&&n!=-1)
{
cout<<solve(n)<<endl;
}
return 0;
}