矩阵快速幂专题（一）

最近闲来无事，准备集中精力刷一波数论与图论。矩阵快速幂是数论里面的重要组成部分，值得我好好学习一下。因为题目比较多，分析也比较多，所以将此专题分成几个部分。做完这一专题，可能会暂时转向图论部分，然后等我组合数学学得差不多了，再回过头来继续做数论题。

矩阵快速幂算法的核心思想是将问题建模转化为数学模型（有一些简单题目是裸的矩阵模型，但是大部分难题就是难在要构造矩阵，用矩阵方法解决问题），推倒递推式，构造计算矩阵，用快速幂的思想求解矩阵A的n次方取mod，从而得到矩阵里面你需要的数据。

矩阵快速幂问题有某些特征，是类似快速幂的形式。首先，从题目中是一定可以得到一个递推式的，而且这个递推式不是简单的An与An-1的关系，而是关系式中带有An-2或之后的项，还有的情况甚至带有常数项和多个数列的组合（加减乘数）；其次，数据量比较大，不可能用打表方法做出来；最后，有些题目可能与矩阵一点没有关系，但是通过规律分析，能够转化为已知递推式求解某一项（取模）的。

这里给出模版：  我觉得我的矩阵模版不是特别好，但是也不是特别麻烦，所以请读者自行选择是否采用

#define maxn 10+5;typedef long long ll;const ll mod = 1e9+7;#define clr(x,y) memset(x,y,sizeof(x))struct matrix {int n;ll maze[maxn][maxn];void init(int n){this->n=n;clr(maze,0);}matrix operator *(const matrix& rhs){matrix ans;ans.init(n);for(int i=0;i<n;i++)for(int j=0;j<n;j++)for(int k=0;k<n;k++)ans.maze[i][j]=(ans.maze[i][j]+maze[i][k]*rhs.maze[k][j])%mod;return ans;}};matrix qlow(matrix a,int n){matrix ans;ans.init(a.n);for(int i=0;i<a.n;i++)ans.maze[i][i]=1;while(n){if(n&1)ans=ans*a;a=a*a;n>>=1;}return ans;}

第一题 hdu-1757 

分析：题目直接给出了递推式（带有f（x-1）到f（x-10）），而且数据量非常大需要取模(这里并不要管a0~a9的值，若不为01那么只需多考虑一步取模)，这是简单的矩阵快速幂裸题。

#include<set>#include<map>#include<cmath>#include<stack>#include<queue>#include<vector>#include<iostream>#include<cstdio>#include<cstring>#include<string>#include<algorithm>#include<cctype>#define maxn 6+5#define clr(x,y) memset(x,y,sizeof(x))using namespace std;const int inf = 0x3f3f3f3f;typedef long long ll;const double pi = acos( -1 );//const ll mod = 1e9+7;int a[10];ll k,m;ll qmul(ll a,ll b){     ll ans=0;    while(b)    {        if(b&1)ans=(ans+a)%m;        a=(a+a)%m;        b>>=1;    }    return ans;}struct matrix{    int n;    ll maze[maxn][maxn];    void init(int n)    {        this->n=n;        clr(maze,0);    }    matrix operator * (matrix & rhs)    {        matrix ans;        ans.init(n);        for(int i=0;i<n;i++)            for(int j=0;j<n;j++)                for(int k=0;k<n;k++)                    ans.maze[i][j]=(ans.maze[i][j]+qmul(maze[i][k],rhs.maze[k][j]))%m;        return ans;    }};matrix qlow(matrix a,ll n){    matrix ans;    ans.init(a.n);    for(int i=0;i<ans.n;i++)ans.maze[i][i]=1;    while(n)    {        if(n&1)ans=ans*a;        a=a*a;        n>>=1;    }    return ans;}int main(){    //freopen("d:\\acm\\in.in","r",stdin);    while(~scanf("%lld %lld",&k,&m))    {        for(int i=0;i<10;i++)            scanf("%d",&a[i]);        if(k<10)        {            printf("%lld\n",k%m);            continue;        }        matrix ans;        ans.init(10);        for(int i=0;i<10;i++)ans.maze[0][i]=i;        matrix ant;        ant.init(10);        for(int i=0;i<9;i++)ant.maze[i+1][i]=1;        for(int i=0;i<10;i++)ant.maze[i][9]=a[9-i];        matrix tmp;        tmp=qlow(ant,k-9);        ans=ans*tmp;        printf("%lld\n",ans.maze[0][9]);    }    return 0;}

第二题 hdu-1575

分析：这道题就是非常直白的矩阵快速幂，将主对角线上的数加起来取模-。-

#include<set>#include<map>#include<cmath>#include<stack>#include<queue>#include<vector>#include<iostream>#include<cstdio>#include<cstring>#include<string>#include<algorithm>#include<cctype>#define maxn 6+5#define clr(x,y) memset(x,y,sizeof(x))using namespace std;const int inf = 0x3f3f3f3f;typedef long long ll;const double pi = acos( -1 );const ll mod = 9973;//1e9+7;ll qmul(ll a,ll b){    ll ans=0;    while(b)    {        if(b&1)ans=(ans+a)%mod;        a=(a+a)%mod;        b>>=1;    }    return ans;}struct matrix{    int n;    ll maze[maxn][maxn];    void init(int n)    {        this->n=n;        clr(maze,0);    }    matrix operator * (matrix& rhs)    {        matrix ans;        ans.init(n);        for(int i=0;i<n;i++)            for(int j=0;j<n;j++)                for(int k=0;k<n;k++)                    ans.maze[i][j]=(ans.maze[i][j]+maze[i][k]*rhs.maze[k][j])%mod;        return ans;    }};matrix qlow(matrix a,int n){    matrix ans;    ans.init(a.n);    for(int i=0;i<ans.n;i++)ans.maze[i][i]=1;    while(n)    {        if(n&1)ans=(ans*a);        a=(a*a);        n>>=1;    }    return ans;}int main(){    //freopen("d:\\acm\\in.in","r",stdin);    int t;    scanf("%d",&t);    while(t--)    {        int n,k;        scanf("%d %d",&n,&k);        matrix ans;        ans.init(n);        for(int i=0;i<n;i++)            for(int j=0;j<n;j++)                scanf("%d",&ans.maze[i][j]);        ans=qlow(ans,k);        ll anw=0;        for(int i=0;i<n;i++)            anw=(anw+ans.maze[i][i])%mod;        printf("%lld\n",anw);    }    return 0;}

第三题 hdu-2604

分析：我这里只提出我的方法，其他网上大神的方法请到别的博客上查找。

对于L较大的情况（L大于2，至于为什么是2，你没有看见fmf和fff都是三个吗？），我考虑对于在E-queue里面的所有合法队列（-。-，其实和队列没有任何关系）。他们都是以mm，fm，mf，ff之一结尾的，我就分别将以他们结尾的并且总长度为n的个数记为mm[n]，fm[n]，mf[n]，ff[n]。

考虑mm[n]后面加m是mm[n+1]的一部分，加f是mf[n+1]的一部分；

fm[n]加f属于O-queue，舍去，加m是mm[n+1]的一部分；

依次类推。。。。。。

得到递推式mm[n+1]=mm[n]+fm[n];    mf[n+1]=mm[n+1];   fm[n+1]=ff[n]+mf[n];    ff[n+1]=mf[n+1];

构造矩阵    

#include<set>#include<map>#include<cmath>#include<stack>#include<queue>#include<vector>#include<iostream>#include<cstdio>#include<cstring>#include<string>#include<algorithm>#include<cctype>#define maxn 0+5#define clr(x,y) memset(x,y,sizeof(x))using namespace std;const int inf = 0x3f3f3f3f;typedef long long ll;const double pi = acos( -1 );const ll mod = 1e9+7;int m;struct matrix{    int n;    ll maze[maxn][maxn];    void init(int n)    {        this->n=n;        clr(maze,0);    }    matrix operator * (const matrix& rhs)    {        matrix ans;        ans.init(n);        for(int i=0;i<n;i++)            for(int j=0;j<n;j++)                for(int k=0;k<n;k++)                    ans.maze[i][j]=(ans.maze[i][j]+maze[i][k]*rhs.maze[k][j])%m;        return ans;    }};matrix qlow(matrix a,int n){    matrix ans;    ans.init(a.n);    for(int i=0;i<ans.n;i++)ans.maze[i][i]=1;    while(n)    {        if(n&1)ans=ans*a;        a=a*a;        n>>=1;    }    return ans;}int main(){    //freopen("d:\\acm\\in.in","r",stdin);    int l;    while(~scanf("%d %d",&l,&m))    {        if(l<=2)        {            if(l==0)puts("0");            else if(l==1)printf("%d\n",2%m);            else printf("%d\n",4%m);            continue;        }        matrix ans;        ans.init(4);        for(int i=0;i<4;i++)ans.maze[0][i]=1;        matrix ant;        ant.init(4);        ant.maze[0][1]=ant.maze[1][3]=ant.maze[2][0]=ant.maze[2][1]=ant.maze[3][2]=ant.maze[3][3]=1;        ant=qlow(ant,l-2);        ans=ans*ant;        int anw=0;        for(int i=0;i<4;i++)anw=(anw+ans.maze[0][i])%m;        printf("%d\n",anw);    }    return 0;}

第四题  hdu-2256

分析：一开始会以为与矩阵无关，但是看见根号还要取模就会发现满满都是套路。-。-虽然根号2与根号3是两个根号，但是还有个2次方，先将2次方乘进去，得到5+2sqrt(6)，发现根号减少到一个，而且对于5+2sqrt(6)的n次方，一定是一个常数项和一个带根号6的项，而且互相之间有递推关系，这就把问题引向矩阵上了。

但是这里不能直接取模，因为Bn还有一部分是sqrt(6)，直接取模显然是错的

这样构造的矩阵和最后要求的数据都有了 -。-

#include<set>#include<map>#include<cmath>#include<stack>#include<queue>#include<vector>#include<iostream>#include<cstdio>#include<cstring>#include<string>#include<algorithm>#include<cctype>#define maxn 0+5#define clr(x,y) memset(x,y,sizeof(x))using namespace std;const int inf = 0x3f3f3f3f;typedef long long ll;const double pi = acos( -1 );const ll mod = 1024;//1e9+7;struct matrix {    int n;    ll maze[2][2];    void init(int n)    {        this->n=n;        clr(maze,0);    }    matrix operator * (const matrix& rhs)    {        matrix ans;        ans.init(n);        for(int i=0;i<n;i++)            for(int j=0;j<n;j++)                for(int k=0;k<n;k++)                    ans.maze[i][j]=(ans.maze[i][j]+maze[i][k]*rhs.maze[k][j])%mod;         return ans;    }};matrix qlow(matrix a,int n){    matrix ans;    ans.init(a.n);    for(int i=0;i<ans.n;i++)ans.maze[i][i]=1;    while(n)    {        if(n&1)ans=ans*a;        a=a*a;        n>>=1;    }    return ans;}int main(){    int t;    scanf("%d",&t);    while(t--)    {        int n;        scanf("%d",&n);        matrix ans;        ans.init(2);        ans.maze[0][0]=1;        matrix ant;        ant.init(2);        ant.maze[0][0]=ant.maze[1][1]=5;        ant.maze[0][1]=2;        ant.maze[1][0]=12;        ant=qlow(ant,n);        ans=ans*ant;        printf("%d\n",(2*ans.maze[0][0]-1)%mod);    }    return 0;}

第五题 codeforces-185A

分析：做完前面的题，这道题反而简单了。不妨考虑当n年后（操作n次），正三角的个数是An，倒三角的个数是Bn。那么继续再操作一次，An个正三角变成了3An个正三角和An个倒三角，Bn个倒三角变成Bn个正三角和3Bn个倒三角。所以有递推关系  An+1=3An+Bn    Bn+1=An+3Bn

构造矩阵   

#include<set>#include<map>#include<cmath>#include<stack>#include<queue>#include<vector>#include<iostream>#include<cstdio>#include<cstring>#include<string>#include<algorithm>#include<cctype>#define maxn 0+5#define clr(x,y) memset(x,y,sizeof(x))using namespace std;const int inf = 0x3f3f3f3f;typedef long long ll;const double pi = acos( -1 );const ll mod = 1e9+7;ll mul(ll a,ll b){ll ans=0;while(b){if(b&1)ans=(ans+a)%mod;a=(a+a)%mod;b>>=1;}return ans;}struct matrix {int n;ll maze[maxn][maxn];void init(int n){this->n=n;clr(maze,0);}matrix operator * (const matrix& rhs){matrix ans;ans.init(n);for(int i=0;i<n;i++)for(int j=0;j<n;j++)for(int k=0;k<n;k++)ans.maze[i][j]=(ans.maze[i][j]+mul(maze[i][k],rhs.maze[k][j]))%mod;return ans;}};matrix qlow(matrix a,ll n){matrix ans;ans.init(a.n);for(int i=0;i<ans.n;i++)ans.maze[i][i]=1;while(n){if(n&1)ans=ans*a;a=a*a;n>>=1;}return ans;}int main(){    //freopen("d:\\acm\\in.in","r",stdin);ll n;scanf("%I64d",&n);matrix ans;ans.init(2);ans.maze[0][0]=1;matrix ant;ant.init(2);ant.maze[0][0]=ant.maze[1][1]=3;ant.maze[1][0]=ant.maze[0][1]=1;ant=qlow(ant,n);ans=ans*ant;printf("%I64d\n",ans.maze[0][0]);    return 0;}

第六题 hdu-2842
分析：题目讲的好像是九连环（好像又不是，英语不好），-。-    可怜我并没有玩过

这道题非常像汉诺塔，但是题目只是求最小步数，而不需要求出具体的步骤。我假设对于长度为i的中国环，我最少需要f(i)步（1<=i<=n）。那么对于长度为n+1的中国环，我们先考虑将最后一个拿下来（没办法，奇葩的规则，只有拿下前n-1个且保留第n个，才能拿下第n+1个）。那么先将前n-1个取下来，需要f(n-1)步；再将第n+1个取下来，需要1步；再把前n-1个放回去，需要f(n-1)步（这里需要解释一下，前n-1个必须放回去，不然第n个无法取下，而且取下需要f(n-1)步，那么放回去也需要f(n-1)，直接相反的步骤再来一遍，不明白的多读几遍题目）；最后，将前n个都取下，需要f(n)步。所以得到递推式 f(n+1)=f(n)+2f(n-1)+1。对于递推式中有常数项，只需要给出一行或列留给常数项。

构造矩阵

#include<set>#include<map>#include<cmath>#include<stack>#include<queue>#include<vector>#include<iostream>#include<cstdio>#include<cstring>#include<string>#include<algorithm>#include<cctype>#define maxn 0+5#define clr(x,y) memset(x,y,sizeof(x))using namespace std;const int inf = 0x3f3f3f3f;typedef long long ll;const double pi = acos( -1 );const ll mod = 200907;//1e9+7;struct matrix {    int n;    ll maze[maxn][maxn];    void init(int n)    {        this->n=n;        clr(maze,0);    }    matrix operator * (const matrix& rhs)    {        matrix ans;        ans.init(n);        for(int i=0;i<n;i++)            for(int j=0;j<n;j++)                for(int k=0;k<n;k++)                    ans.maze[i][j]=(ans.maze[i][j]+maze[i][k]*rhs.maze[k][j])%mod;        return ans;    }};matrix qlow(matrix a,int n){    matrix ans;    ans.init(a.n);    for(int i=0;i<ans.n;i++)ans.maze[i][i]=1;    while(n)    {        if(n&1)ans=ans*a;        a=a*a;        n>>=1;    }    return ans;}int main(){    //freopen("d:\\acm\\in.in","r",stdin);    int n;    matrix ant;    ant.init(3);    ant.maze[0][0]=ant.maze[0][1]=ant.maze[2][0]=ant.maze[2][2]=1;    ant.maze[1][0]=2;    while(scanf("%d",&n),n)    {        if(n<=2)        {            if(n==1)printf("%d\n",1);            else printf("%d\n",2);            continue;        }        matrix ans;        ans.init(3);        ans.maze[0][0]=2;        ans.maze[0][1]=ans.maze[0][2]=1;        ans=ans*qlow(ant,n-2);        printf("%I64d\n",ans.maze[0][0]);    }    return 0;}

第七题 hdu-2276

分析：对于01串，如果前一位的数位为1，那么这一位就翻转（0变1，1变0）；前一位为0，那么这一位就不变。这就像前一位与这一位相加，并且mod2，如果前一位是0，那么这一位加上前一位再mod2数值不会改变；如果前一位是1，那么这一位就相当于+1mod2，正好和翻转的定义相同。这样我们就可以构造出矩阵

#include<set>#include<map>#include<cmath>#include<stack>#include<queue>#include<vector>#include<iostream>#include<cstdio>#include<cstring>#include<string>#include<algorithm>#include<cctype>#define maxn 100+5#define clr(x,y) memset(x,y,sizeof(x))using namespace std;const int inf = 0x3f3f3f3f;typedef long long ll;const double pi = acos( -1 );const ll mod = 1e9+7;char str[maxn];struct matrix {    int n;    ll maze[maxn][maxn];    void init(int n)    {        this->n=n;        clr(maze,0);    }    matrix operator * (const matrix & rhs)    {        matrix ans;        ans.init(n);        for(int i=0;i<n;i++)            for(int j=0;j<n;j++)                for(int k=0;k<n;k++)                    ans.maze[i][j]=(ans.maze[i][j]+maze[i][k]*rhs.maze[k][j])%2;        return ans;    }};matrix qlow(matrix a,int n){    matrix ans;    ans.init(a.n);    for(int i=0;i<ans.n;i++)ans.maze[i][i]=1;    while(n)    {        if(n&1)ans=ans*a;        a=a*a;        n>>=1;    }    return ans;}int main(){    //freopen("d:\\acm\\in.in","r",stdin);    int m;    while(~scanf("%d",&m))    {        scanf("%s",str);        int n=strlen(str);        matrix ans;        ans.init(n);        for(int i=0;i<n;i++)            if(str[i]=='1')ans.maze[0][i]=1;        matrix ant;        ant.init(n);        for(int i=0;i<n;i++)        {            ant.maze[i][i]=1;            if(i==0)ant.maze[n-1][0]=1;            else ant.maze[i-1][i]=1;        }        ant=qlow(ant,m);        ans=ans*ant;        for(int i=0;i<n;i++)            printf("%d",ans.maze[0][i]);        puts("");    }    return 0;}

第八题 hdu-2254

分析：学过离散数学的都知道，如果用邻接矩阵来存图的话，矩阵A的n次方里面的（A^n）[i][j]表示从i到j且路径长度为n的不同路径条数。这一定理正好适用这一题目。我用邻接矩阵来建图（单向边，一条边对应一个++操作）得到矩阵A，题目要求的就是｛A^(t1)｝[v1][v2]+｛A^(t1+1)｝[v1][v2]+｛A^(t1+2)｝[v1][v2]+......+｛A^(t2)｝[v1][v2]的和取模。但是如果一个一个算的话特别麻烦，还有超时的可能，所以我们进一步改进思路。

对于，我们不妨设Z(n)=A^1+A^2+...+A^n

那么可以构造递推关系 A^(n+1)=A^n * A      Z(n+1)=Z(n)+ A^(n+1) =Z(n)+A^n * A 

构造矩阵

但是这道题目略坑，有几个注意点：

1、城市的编号不是从0到n-1，而是随便的一个数字，需要离散化否则不能存相关信息

2、城市数不超过30，也就是说我的方法开矩阵不超过60，但是我残念的一开始以为最多可能有20000个不同城市    血崩！

3、图中可能有重边，所以别用=1，要用++操作

4、询问中v1，v2可能在前面的城市编号集中没有出现，那么此时答案为0

5、t1可能比t2大，这种情况你就交换下t1，t2

这道题要做对也是真心不容易！

#include<set>#include<map>#include<cmath>#include<stack>#include<queue>#include<vector>#include<iostream>#include<cstdio>#include<cstring>#include<string>#include<algorithm>#include<cctype>#define maxn 60+5#define clr(x,y) memset(x,y,sizeof(x))using namespace std;const int inf = 0x3f3f3f3f;typedef long long ll;const double pi = acos( -1 );const ll mod = 2008;//1e9+7;ll gra[10000+5][2];ll num[20000+5];int sum;int bs(ll k){    int m,l=0,r=sum-1;    while(l<=r)    {        m=(l+r)>>1;        if(k==num[m])return m;        else if(k>num[m])l=m+1;        else r=m-1;    }    return -1;}struct matrix {    int n;    ll maze[maxn][maxn];    void init(int n)    {        this->n=n;        clr(maze,0);    }    matrix operator * (const matrix & rhs)    {        matrix ans;        ans.init(n);        for(int i=0;i<n;i++)            for(int j=0;j<n;j++)                for(int k=0;k<n;k++)                    ans.maze[i][j]=(ans.maze[i][j]+maze[i][k]*rhs.maze[k][j])%mod;        return ans;    }};matrix qlow(matrix a,int n){    matrix ans;    ans.init(a.n);    for(int i=0;i<a.n;i++)ans.maze[i][i]=1;    while(n)    {        if(n&1)ans=ans*a;        a=a*a;        n>>=1;    }    return ans;}int main(){    //freopen("d:\\acm\\in.in","r",stdin);    int n;    while(~scanf("%d",&n))    {        sum=n*2;        for(int i=0;i<n;i++)        {            scanf("%lld %lld",&gra[i][0],&gra[i][1]);            num[i<<1]=gra[i][0];            num[i<<1|1]=gra[i][1];        }        sort(num,num+sum);        sum=unique(num,num+sum)-num;        matrix ans;        ans.init(sum<<1);        matrix ant;        ant.init(sum<<1);        int a,b;        for(int i=0;i<n;i++)        {            a=bs(gra[i][0]);            b=bs(gra[i][1]);            ans.maze[a][b]++;            ans.maze[a][b+sum]++;        }        for(int i=0;i<sum;i++)            ans.maze[i+sum][i+sum]=1;        int q;        scanf("%d",&q);        while(q--)        {            int t1,t2;            ll v1,v2;            scanf("%lld %lld %d %d",&v1,&v2,&t1,&t2);            v1=bs(v1);            v2=bs(v2);            if(v1==-1||v2==-1)            {                puts("0");                continue;            }            if(t1>t2)swap(t1,t2);            t1--;            int d,d1,d2;            if(t1<=0)d1=0;            else             {                ant=qlow(ans,t1);                d1=ant.maze[v1][v2+sum];            }            if(t2<=0)d2=0;            else             {                ant=qlow(ans,t2);                d2=ant.maze[v1][v2+sum];            }            d=d2-d1;            while(d<0)d+=mod;            printf("%d\n",d);        }    }    return 0;}

第九题 hdu-3117

分析：其实这道题不应该放在矩阵快速幂专题的，因为这道题关于矩阵的部分很简单-。-，难的是题目的数学部分-。-

求解斐波那契数列的值，对于数值位不超过8位的直接输出数值；但是对于大于8位的，输出值的前四位和后四位（格式自己看-。-）。后四位当然非常简单了，f(n)=f(n-1)+f(n-2) ,mod 10000,我就不再赘述了

但是对于前四位还是挺麻烦的，这里要用到对数的优良性质。logA(b^c)=c*logA(c)   logA(b*c)=logA(b)+logA(c)

对于一个大数（不如就选10234432）,怎么知道它的前四位呢

log10(10234432)=log10(1.0234432*10^7)=log10(1.0234432)+7

7是log10(10234432)的整数部分，log10(1.0234432)是log10(10234432)的小数部分

可以明显看出pow(10,log10(1.0234432))=1.0234432   从而可以简单地得到了前几个位数

有人会问为什么要这么麻烦，不如以字符串的形式读入然后输出前四位不就好了吗 ？

对！没错，对于确定的大数，确实是可以的，但是对b^c来说却不能使用，但是取对数的方法却可以

对于这道题  我们考虑斐波那契数列的通项（有人问，比赛的时候你也知道通项？我的回答是可以，给我类似斐波那契的递推式我就可以求出其通项，这个在数学竞赛中有学到，请读者自行百度）

F(n)=(1/sqrt(5))*{[(1+sqrt(5))/2]^n-[(1-sqrt(5))/2]^n}=(1/sqrt(5))*{[(1+sqrt(5))/2]^n}*{1-[(1-sqrt(5))/(1+sqrt(5))]^n}

设p=(1-sqrt(5))/2

s=log10(F(n))=-1/2*log10(5)+n*log10(p)+Δ（Δ非常小，可以忽略不计，Δ小于log10(1)）

s去掉它的整数部分，然后(int) [ pow（10，s-（int）s）*1000 ] 就得到了结果

如果想看更加详细的分析请点大神的博客 http://blog.sina.com.cn/s/blog_9bf748f301019q3t.html

我只是简单搬运而且=、=

#include<set>#include<map>#include<cmath>#include<stack>#include<queue>#include<vector>#include<iostream>#include<cstdio>#include<cstring>#include<string>#include<algorithm>#include<cctype>#define maxn 60+5#define clr(x,y) memset(x,y,sizeof(x))using namespace std;const int inf = 0x3f3f3f3f;typedef long long ll;const double pi = acos( -1 );const ll mod = 10000;//1e9+7;int fa[60];struct matrix {    int n;    ll maze[5][5];    void init(int n)    {        this->n=n;        clr(maze,0);    }    matrix operator * (const matrix & rhs)    {        matrix ans;        ans.init(n);        for(int i=0;i<n;i++)            for(int j=0;j<n;j++)                for(int k=0;k<n;k++)                    ans.maze[i][j]=(ans.maze[i][j]+maze[i][k]*rhs.maze[k][j])%mod;        return ans;    }};matrix qlow(matrix a,int n){    matrix ans;    ans.init(a.n);    for(int i=0;i<a.n;i++)ans.maze[i][i]=1;    while(n)    {        if(n&1)ans=ans*a;        a=a*a;        n>>=1;    }    return ans;}int main(){    //freopen("d:\\acm\\in.in","r",stdin);    fa[0]=0;    fa[1]=1;    for(int i=2;i<40;i++)        fa[i]=fa[i-1]+fa[i-2];    int n;    while(~scanf("%d",&n))    {        if(n<40)        {            printf("%d\n",fa[n]);            continue;        }        double p=(sqrt(5.0)+1.0)/2.0;        double s=1.0*n*log10(p)-log10(5.0)/2;        s=s-(int)s;        int anw=(int)(pow(10.0,s)*1000);        printf("%d...",anw);        matrix ant;        ant.init(2);        ant.maze[0][0]=ant.maze[0][1]=ant.maze[1][0]=1;        ant=qlow(ant,n-1);        printf("%04d\n",ant.maze[0][0]);    }    return 0;}

第十题 hdu-4686

分析：这是一道显然的矩阵题（刷完前面那么多道，如果还没有一眼看出来的话，建议停下来重新理解算法-。-）

感觉难度不是很大，但是这道题有多个递推关系，而且还涉及到常数项和两个不相关的递推关系乘积，稍微有点麻烦。但是冷静下来，对每一个递推关系进行分析和对你所要求的项进行分解。

由AoD(n)=a(0)b(0)+a(1)b(1)+...+a(n-1)b(n-1)   可得    AoD(n+1)=AoD(n)+a(n)b(n)

那么a(n)b(n)怎么求呢   现在只知道  a(n)=a(n-1)*ax+ay    b(n)=b(n-1)*bx+by

那么可不可以利用已经得到的a(n-1)和b(n-1)和a(n-1)b(n-1)   

可以得到 a(n)b(n)=ax*bx*a(n-1)*b(n-1)+ax*by*a(n-1)+ay*bx*b(n-1)+ay*by

那么问题又来了  a(n)和 b(n)又怎么单独求呢     这个题目已经给出  a(n)=a(n-1)*ax+ay    b(n)=b(n-1)*bx+by

分析清楚了所有的递推关系后，可以来构造矩阵了

注意点：a0，ax，ay，b0，bx，by都可能大于1e7+9，所以需要先取一次模

#include<set>#include<map>#include<cmath>#include<stack>#include<queue>#include<vector>#include<iostream>#include<cstdio>#include<cstring>#include<string>#include<algorithm>#include<cctype>#define maxn 1+5#define clr(x,y) memset(x,y,sizeof(x))using namespace std;const int inf = 0x3f3f3f3f;typedef long long ll;const double pi = acos( -1 );const ll mod = 1e9+7;ll mul(ll a,ll b){    ll ans=0;    while(b)    {        if(b&1)ans=(ans+a)%mod;        a=(a+a)%mod;        b>>=1;    }    return ans;}struct matrix {    int n;    ll maze[maxn][maxn];    void init(int n)    {        this->n=n;        clr(maze,0);    }    matrix operator * (const matrix & rhs)    {        matrix ans;        ans.init(n);        for(int i=0;i<n;i++)            for(int j=0;j<n;j++)                for(int k=0;k<n;k++)                    ans.maze[i][j]=(ans.maze[i][j]+maze[i][k]*rhs.maze[k][j])%mod;        return ans;    }};matrix qlow(matrix a,ll n){    matrix ans;    ans.init(a.n);    for(int i=0;i<a.n;i++)ans.maze[i][i]=1;    while(n)    {        if(n&1)ans=ans*a;        a=a*a;        n>>=1;    }    return ans;}int main(){    //freopen("d:\\acm\\in.in","r",stdin);    ll n;    while(~scanf("%lld",&n))    {        ll a0,b0,ax,bx,ay,by;        scanf("%lld %lld %lld %lld %lld %lld",&a0,&ax,&ay,&b0,&bx,&by);        a0%=mod;        b0%=mod;        ax%=mod;        bx%=mod;        ay%=mod;        by%=mod;        if(n==0)        {            puts("0");            continue;        }        matrix ans;        ans.init(5);        ans.maze[0][0]=1;        ans.maze[0][1]=a0;        ans.maze[0][2]=b0;        ans.maze[0][3]=a0*b0%mod;        matrix ant;        ant.init(5);        ant.maze[0][0]=1;        ant.maze[0][1]=ay;        ant.maze[1][1]=ax;        ant.maze[0][2]=by;        ant.maze[2][2]=bx;        ant.maze[0][3]=ay*by%mod;        ant.maze[1][3]=ax*by%mod;        ant.maze[2][3]=ay*bx%mod;        ant.maze[3][3]=ax*bx%mod;        ant.maze[3][4]=ant.maze[4][4]=1;        ant=qlow(ant,n);        ans=ans*ant;        printf("%lld\n",ans.maze[0][4]);    }    return 0;}

这里先小结一下，矩阵快速幂类型的题目非常像科学归纳法（第一科学归纳法、第二科学归纳法还有其他类型的），主要是求得当前项与前几项或者常数项的递推关系，从而构造矩阵求解。