矩阵（matrix）应用大总结（二）HDU1588+HDU 5950+HDU 5015

        聪神选的题系列……

        这一系列可以说是数学系列咯，正好补一补。

        这几道题目都是矩阵的题目，本来呢觉得自己之前也学过见过挺多矩阵的，还自己做过总结，做这些题目应该没什么大问题，但是还是感觉有很多新东西。

        首先，我们说说 HDU 1588 吧，Gauss Fibonacci数列，就是求f(g(i))的和，其中f(i)表示fibonacci数列，g(i)=k*i+b(0<=i<n)。换句话说就是求f(b)+f(k+b)+f(2k+b)+..+f(nk+b)。然后n、k、b的范围都是1e9。根据之前总结的，求fibonacci数列的某一项以及它的前n项和很简单，但是这样求离散的几项的和怎么求呢？当然不要那么死板嘛，当初我们求前n项和利用的是转移矩阵等比的性质，那么这里可以吗？我们发现f(k+b)=f(b)*A^k，f(2k+b)=f(b)*A^2k，...，f((n-1)k+b)=f(b)*A^(n-1)k，其中A为递推的转移矩阵。这样，相当于求f(b)*(1+A^k+A^2k+...+A^(n-1)k)，也是一个等比矩阵求和，用同样的方法解决即可，但是在要在原模板上加以修改，具体见代码：

#include<bits/stdc++.h>#define LL long long#define N 3using namespace std;LL n,m,k,b,mod;struct matrix{LL a[N][N];void init(){for(int i=1;i<=n;i++)for(int j=1;j<=n;j++)a[i][j]=0;}} x,basic;matrix matmulti(matrix x,matrix y){matrix ans; ans.init();for(int i=1;i<=n;i++)for(int j=1;j<=n;j++){for(int k=1;k<=n;k++)        ans.a[i][j]+=(x.a[i][k]*y.a[k][j])%mod;    ans.a[i][j]%=mod;}return ans;}matrix matadd(matrix x,matrix y){for(int i=1;i<=n;i++)for(int j=1;j<=n;j++)    x.a[i][j]=(x.a[i][j]+y.a[i][j])%mod;return x;}matrix matpow(matrix x,LL y){    if (y==0) return basic;while ((y&1)==0){y>>=1;x=matmulti(x,x);}matrix ans=x; y>>=1;while (y!=0){x=matmulti(x,x);if ((y&1)!=0) ans=matmulti(ans,x);y>>=1;}return ans;}matrix matmultisum(matrix x,LL y){if (y==1) return matpow(x,b);//所有的额项都要乘上一个f(b)故返回的是x^bmatrix ans;ans=matmultisum(x,y>>1);matrix z=matpow(x,(y>>1)*k);//等比矩阵的公比是x^kans=matadd(ans,matmulti(ans,z));if (y&1) ans=matadd(ans,matpow(x,(y-1)*k+b));return ans;}int main(){    n=2;    basic.a[1][1]=basic.a[2][2]=1;    x.a[1][1]=x.a[1][2]=x.a[2][1]=1;    while (~scanf("%I64d%I64d%I64d%I64d",&k,&b,&m,&mod))    {        matrix ans=matmultisum(x,m);        printf("%I64d\n",ans.a[1][2]);    }    return 0;}

        然后我们说回矩阵的递推。我在第一次总结的时候给出了一种线性递推矩阵的基本构造方法（其实是叉姐和许多大神给出的），现在我要打破原来的方法，并同时把递推从线性推到非线性。首先，说说本质，矩阵递推本质是使一个向量乘以转移矩阵而得到下一个向量。通常这个向量里放的就是所有与递推有关的项。以fibonacci数列为例，初始向量为(an-1,an-2)^-1，一次转移以后需要的到的向量是(an,an-1)^-1。现在推转移矩阵A，根据矩阵乘法的性质，A的第一行乘以向量得到an，A的第二行乘以向量得到an-1，即A(1,1)*an-1+A(1,2)*an-2=an 和 A(2,1)*an-1+A(2,2)*an-2=an-1。根据递推关系，可以退出A(1,1)=A(1,2)=A(2,1)=1，A(2,2)=0。这样子我们就把矩阵给构造出来了，对于任意的递推关系都是如此（有空可以用这种方式退一下Lost In WHU那题的转移矩阵）。差点忘说了，这样子的话我们需要的结果是转移矩阵的幂的矩阵的第一行乘以初始值，还是一fibonacci数列为例，f(n)=A^(n-1)(1,1)*f(1)+A^(n-1)(1,2)*f(0)。按照这个来说我上一题的写法还不规范……

        接着，用一个非线性递推的例子巩固一下这种构造方式。HDU 5950 。题意不说了，就是f(n)=2*f(n-2)+f(n-1)+n^4，然后求第n项。你没有看错，就是n^4，这可是去年沈阳站regional的题目。对于这种非线性的东西，我们一定要学会分解，将非线性转化为线性的，或者说转换成也有递推关系的东西。注意到n^4=(n-1+1)^4=(n-1)^4+C(4,1)(n-1)^3+C(4,2)(n-1)^2+C(4,3)(n-1)+1。很漂亮的，n^4就和四个项有了递推关系，而且这四个项刚好变量都为(n-1)。于是初始向量为(f(n-1),f(n-2),(n-1)^4,(n-1)^3,(n-1)^2,n-1,1)^-1，一次转移后需要得到的向量是(f(n),f(n-1),n^4,n^3,n^2,n,1)。n^4的转移方法刚刚说了，n^3、n^2的转移方法也和n^4类似，也是分解。按照矩阵乘法的性质，我们可以得到下图的矩阵（按照我的习惯，其中图中所有向量和矩阵的第一行和第二行需要交换一下位置）：

                             

        可以说这幅图很好的诠释了我所说的矩阵构造法，就是先写出初始向量，和下一个向量，把所以存在在递推式中的东西都放到向量里面，之后根据矩阵乘法的定义，一个个的求出转移矩阵每个地方应该填什么东西。下面给出具体代码：

#include<bits/stdc++.h>#define LL long long#define mod 2147493647#define N 8using namespace std;struct matrix{LL a[N][N];void init(){for(int i=1;i<N;i++)for(int j=1;j<N;j++)a[i][j]=0;}} x,basic;int n=7;LL y;matrix matmulti(matrix x,matrix y){matrix ans; ans.init();for(int i=1;i<=n;i++)for(int j=1;j<=n;j++){for(int k=1;k<=n;k++)        ans.a[i][j]+=(x.a[i][k]*y.a[k][j])%mod;    ans.a[i][j]%=mod;}return ans;}matrix matpow(matrix x,LL y){if (y==0) return basic;elsewhile ((y&1)==0){y>>=1;x=matmulti(x,x);}matrix ans=x; y>>=1;while (y!=0){x=matmulti(x,x);if ((y&1)!=0) ans=matmulti(ans,x);y>>=1;}return ans;}int main(){    int T_T;    x.a[1][2]=2;    scanf("%d",&T_T);    x.a[1][1]=x.a[2][1]=1;    x.a[1][7]=x.a[1][3]=1;    for(int i=3;i<=7;i++)        x.a[i][i]=x.a[i][7]=1;    x.a[3][4]=x.a[1][4]=4;    x.a[3][5]=x.a[1][5]=6;    x.a[3][6]=x.a[1][6]=4;    x.a[4][5]=3; x.a[4][6]=3; x.a[5][6]=2;    while (T_T--)    {        LL m,a,b;        scanf("%I64d%I64d%I64d",&m,&a,&b);        if (m==1) printf("%I64d\n",a);        else if (m==2) printf("%I64d\n",b);        else        {            matrix ans=matpow(x,m-2);            printf("%I64d\n",(ans.a[1][1]*b%mod+ans.a[1][2]*a%mod+ans.a[1][3]*16+ans.a[1][4]*8+ans.a[1][5]*4+ans.a[1][6]*2+ans.a[1][7])%mod);//最后的数值是转移矩阵的幂的第一行乘以初始值向量        }    }    return 0;}

        矩阵以及矩阵构造的巧妙不仅这些。再看看 HDU 5015 。第一行全是233……，下面几行给出a(i,0)，然后递推式a(i,j)=a(i-1,j)+a(i,j-1)。对，很像组合数。说实话，看到了很懵逼，二维递推，如果可以的话不久可以用来求组合数了吗，为什么之前学的时候没听过……其实是有限制的，这里n<10，关于求组合数的问题等等会说明一下。二维递推的话还是一样的找关系，由于实在不可能用二维的递推关系，所以我们以一列为一个单位(n<10，可承受)。初始向量(0,a(1,0),a(2,0),...,a(n,0))，转移一次之后的向量(233,a(1,0)+233,a(2,0)+a(1,0)+233,...,a(n,0)+a(n-1,0)+...+a(1,0)+233)。看起来很复杂，于是我们把首项的0变成23，然后由于第一行的数字有a(0,m)=a(0,m-1)*10+3的关系，我们把10和3也写到向量中，再交换一下次序，用m代替具体数字则有：初始向量(a(m-1,n-1),a(m-2,n-1),..,a(0,n-1),10,3)，转移一次之后的向量(a(m,n-1),a(m-1,n-1),...,a(0,n-1),10,3)。这里我把刚刚说的递推关系再说清楚点（以n=4为例）：

第一列：0(23)              第二列：23*10+3
              a1                                 23*10+3+a1
              a2                                 23*10+3+a1+a2
              a3                                 23*10+3+a1+a2+a3
              a4                                 23*10+3+a1+a2+a3+a4

        这样应该很清楚明了了，我只不过是按照我个人的习惯换了一下顺序。这样的话就可以得到如下的矩阵（还是以n=4为例）：

                                                            

        一定要注意，我这个顺序是换了个位置的，因为我个人的习惯，最后的结果一定是转移矩阵的幂的第一行乘以初始首向量。具体见代码：

#include<bits/stdc++.h>#define LL long long#define mod 10000007#define N 13using namespace std;struct matrix{LL a[N][N];void init(){for(int i=1;i<N;i++)for(int j=1;j<N;j++)a[i][j]=0;}} x,basic;int n,m,a[11];matrix matmulti(matrix x,matrix y){matrix ans; ans.init();for(int i=1;i<=n;i++)for(int j=1;j<=n;j++){for(int k=1;k<=n;k++)        ans.a[i][j]+=x.a[i][k]*y.a[k][j];    ans.a[i][j]%=mod;}return ans;}matrix matpow(matrix x,LL y){if (y==0) return basic;elsewhile ((y&1)==0){y>>=1;x=matmulti(x,x);}matrix ans=x; y>>=1;while (y!=0){x=matmulti(x,x);if ((y&1)!=0) ans=matmulti(ans,x);y>>=1;}return ans;}int main(){    for(int i=1;i<N;i++)        basic.a[i][i]=1;    while (~scanf("%d%d",&n,&m))    {        for(int i=1;i<=n;i++)            scanf("%d",&a[i]);        if (m==0)        {            printf("%d\n",a[n]);            continue;        } n+=2;        for(int i=1;i<=n;i++)            for(int j=i;j<=n;j++)                x.a[i][j]=1;//构造转移矩阵        for(int i=1;i<n;i++) x.a[i][n-1]=10;        matrix ans=matpow(x,m);        int res=0;        for(int i=1;i<=n-2;i++)            res+=ans.a[1][i]*a[n-i-1]%mod;//结果一定是转移矩阵的幂的第一行乘以初始首向量        printf("%d\n",(res+ans.a[1][n-1]*23+ans.a[1][n]*3)%mod);    }    return 0;}