矩阵十题（7）

Matrix67原文：经典题目7 VOJ1067

我们可以用上面的方法二分求出任何一个线性递推式的第n项，其对应矩阵的构造方法为：在右上角的(n-1)*(n-1)的小矩阵中的主对角线上填1，矩阵第n行填对应的系数，其它地方都填0。例如，我们可以用下面的矩阵乘法来二分计算f(n) = 4f(n-1) - 3f(n-2) + 2f(n-4)的第k项：

利用矩阵乘法求解线性递推关系的题目我能编出一卡车来。这里给出的例题是系数全为1的情况。

下面是题解

vijos 1067 Warcraft III 守望者的烦恼

题目链接：https://vijos.org/p/1067

题目大意：有n个监狱排成一列，每次最多可以往前走k个监狱，问走到第n个监狱有多少种方案，结果mod 7777777

1<=k<=10, 1<=n<=2^31-1

这道题显然是一道DP，但是数据范围太大了，所以我们要去优化整个DP。

设f[n]为监狱长度为n时的方案数，令f[0]=1

我们从起点出发，每次最多可以走k步，求我们到达点n的方案数f[n]，可以分成刚开始走1步，2步，3步....k步的情况，如果刚开始走的是1步，则问题变成走了1步之后，我们从起点出发到达点n-1的方案数f[n-1]，

如果刚开始走的是2步，则问题变成走了2步之后，我们从起点出发到达点n-2的方案数f[n-2]，

以此类推.....

如果刚开始走的是k步，则问题变成走了k步之后，我们从起点出发到达点n-k的方案数f[n-k]，

我们可以发现，这就变成了一个递推的问题了，然后只要我们把走1步，走2步.....走k步后对应的方案数加起来就变成我们要求的f[n]了，

即状态转移方程为：f[n]=f[n-1]+f[n-2]+....f[n-k]。

f[0]=1

当k=1, f[1]=1;

f[2]=f[1]=1;

f[3]=f[2]=1;

f[4]=f[3]=1;

当k=2, f[1]=1;

f[2]=f[1]+f[0]=2;

f[3]=f[2]+f[1]=3;

f[4]=f[3]+f[2]=5;

当k=3, f[1]=1;

f[2]=f[1]+f[0]=2;

f[3]=f[2]+f[1]+f[0]=4;

f[4]=f[3]+f[2]+f[1]=7;

我们手算验证一下发现符合递推式，故推出来的状态转移方程是正确的。

其实任何一个递推式都可以用矩阵乘法来进行优化，我们可以用二分求出任何一个线性递推式的第n项，

其对应矩阵的构造方法为：（和Matrix67的方法有点小不同，此处是Matrix67方法构造的矩阵的转置）

在左下角的(n-1)*(n-1)的小矩阵中的主对角线上填1，矩阵第n列填对应的系数，其它地方都填0。

例如，我们可以用下面的矩阵乘法来二分计算f(n) = 4f(n-1) - 3f(n-2) + 2f(n-4)的第k项：

[ [ 0 1 0 0] [f(k-4)] [f(k-3)]

[ 0 0 1 0] * [f(k-3)] = [f(k-2)]

[ 0 0 0 1] [f(k-2)] [f(k-1)]

[ 2 0 -3 4] ] [f(k-1)] [f( k )]

观察一下DP转移方程：F[i]=F[i-1]+F[i-2]+...+F[i-k];注意到这是一个形如f[i]=a1*f[i-1]+a2*f[i-2]+...+ak[i-k];的方程，立刻就反应到可以使用矩阵乘法进行优化。

因此构造一个k*k的矩阵A然后求n-k次幂A^(n-k)，然后用递推式构造前k项，然后用之与前面求出的矩阵A^(n-k)相乘得到f[n]即可。

我的构造方法，以k=2为例

[f(1),f(2)] * [ [ 0 1 ] ^(n-k) = [f(n-1),f(n)]

[ 1 1 ] ]

代码如下：

 1 #include<stdio.h>
 2 #include<string.h>
 3 #define N 11
 4 #define M 7777777
 5 struct Matrix
 6 {
 7     __int64 a[N][N];
 8 }res,tmp,origin,A,ans;
 9 int n,m,f[N];
10 Matrix mul(Matrix x,Matrix y)
11 {
12     int i,j,k;
13     memset(tmp.a,0,sizeof(tmp.a));
14     for(i=1;i<=n;i++)
15         for(j=1;j<=n;j++)
16             for(k=1;k<=n;k++)
17             {
18                 tmp.a[i][j]+=(x.a[i][k]*y.a[k][j])%M;
19                 tmp.a[i][j]%=M;
20             }
21     return tmp;
22 }
23 void quickpow(int k)  //矩阵快速幂
24 {
25     int i;
26     memset(res.a,0,sizeof(res.a));
27     for(i=1;i<=n;i++)
28         res.a[i][i]=1;
29     while(k)
30     {
31         if(k&1)
32             res=mul(res,A);
33         A=mul(A,A);
34         k>>=1;
35     }
36 }
37 int main()
38 {
39     int k,i,j;
40     while(scanf("%d%d",&k,&m)!=EOF)
41     {
42         memset(f,0,sizeof(f));
43         f[0]=1;
44         for(i=1;i<=k;i++)      //构造前k项
45             for(j=0;j<i;j++)
46                 f[i]+=f[j];
47         memset(A.a,0,sizeof(A.a));
48         for(i=2;i<=k;i++)      //构造矩阵A
49             A.a[i][i-1]=1;
50         for(i=1;i<=k;i++)
51             A.a[i][k]=1;
52         n=k;
53         quickpow(m-k);   //A^(n-k)
54         memset(origin.a,0,sizeof(origin.a));
55         for(i=1;i<=k;i++)     //前k项的矩阵[f(1),f(2),f(3)....f(k)]
56             origin.a[1][i]=f[i];
57         ans=mul(origin,res);      //[f(1),f(2),f(3)....f(k)] * A^(n-k)
58         printf("%d\n",ans.a[1][n]%M);
59     }
60     return 0;
61 }

View Code

hdu 1757 A Simple Math Problem

题目链接：http://acm.hdu.edu.cn/showproblem.php?pid=1757

If x < 10 ，则 f(x) = x.
If x >= 10 ，则 f(x) = a0 * f(x-1) + a1 * f(x-2) + a2 * f(x-3) + …… + a9 * f(x-10);

给出k，m和a0~a9，求f(k)%m， k<2*10^9 , m < 10^5

这是一个递推式，故可以用矩阵乘法来求

根据上面讲的方法构建矩阵A,然后二分求A^(k-10),然后用前k项的矩阵与之相乘即可得到f[n]

代码如下：

 1 #include<stdio.h>
 2 #include<string.h>
 3 #define N 15
 4 struct Matrix
 5 {
 6     __int64 a[N][N];
 7 }ans,tmp,res,origin,A;
 8 int n,m,f[N],a[N];
 9 Matrix mul(Matrix x,Matrix y)
10 {
11     int i,j,k;
12     memset(tmp.a,0,sizeof(tmp.a));
13     for(i=1;i<=n;i++)
14         for(j=1;j<=n;j++)
15             for(k=1;k<=n;k++)
16             {
17                 tmp.a[i][j]+=(x.a[i][k]*y.a[k][j])%m;
18                 tmp.a[i][j]%=m;
19             }
20     return tmp;
21 }
22 Matrix quickpow(int k)  //矩阵快速幂
23 {
24     int i;
25     memset(res.a,0,sizeof(res.a));
26     for(i=1;i<=n;i++)
27         res.a[i][i]=1;
28     while(k)
29     {
30         if(k&1)
31             res=mul(res,A);
32         A=mul(A,A);
33         k>>=1;
34     }
35     return res;
36 }
37 int main()
38 {
39     int k,i;
40     while(scanf("%d%d",&k,&m)!=EOF)
41     {
42         n=10;
43         for(i=0;i<n;i++)
44             scanf("%d",&a[i]);
45         f[n]=0;
46         for(i=0;i<n;i++)   //求前10项
47             f[i]=i;
48         for(i=0;i<n;i++)
49             f[n]+=a[i]*f[n-1-i];
50         if(k<=n)
51             printf("%d\n",f[k]);
52         else
53         {
54             memset(A.a,0,sizeof(A.a));
55             for(i=2;i<=n;i++)    //构造矩阵A
56                 A.a[i][i-1]=1;
57             for(i=1;i<=n;i++)
58                 A.a[i][n]=a[n-i];
59             res=quickpow(k-10);   //A^(k-10)
60             memset(origin.a,0,sizeof(origin.a));
61             for(i=1;i<=n;i++)   //前10项的矩阵[f(1),f(2),f(3)....f(10)]
62                 origin.a[1][i]=f[i];
63             ans=mul(origin,res);    //[f(1),f(2),f(3)....f(10)] * A^(n-10)
64             printf("%d\n",ans.a[1][n]%m);
65         }
66     }
67     return 0;
68 }

View Code