利用Cayley-Hamilton theorem 优化矩阵线性递推
Posted lzw4896s
tags:
篇首语：本文由小常识网(cha138.com)小编为大家整理，主要介绍了利用Cayley-Hamilton theorem 优化矩阵线性递推相关的知识，希望对你有一定的参考价值。
平时有关线性递推的题,很多都可以利用矩阵乘法来解决。 时间复杂度一般是O(K3logn)因此对矩阵的规模限制比较大。
下面介绍一种利用利用Cayley-Hamilton theorem加速矩阵乘法的方法。
Cayley-Hamilton theorem：
记矩阵A的特征多项式为f(x)。 则有f(A)=0.
证明可以看 维基百科 https://en.wikipedia.org/wiki/Cayley–Hamilton_theorem#A_direct_algebraic_proof
另外我在高等代数的课本上找到了 证明(和维基百科里的第一种证明方法是一样的)
下面介绍几个 利用可以这个定理解决的题目：
显然可以用矩阵乘法来做。下面讲一下怎么利用Cayley-Hamilton theorem 来优化： 详细的论述可以参考这篇。
设M为K阶矩阵,主要思想就是将$M^n$表示为 $b_0 M^0\\ +\\ b_1 M^1\\ +\\ \\cdots\\ b_{K-1}M^{K-1}$这样的形式.
根据Cayley-Hamilton theorem $M^K\\ =\\ a_0 M^0\\ +\\ a_1 M^1\\ +\\ \\cdots\\ a_{K-1}M^{K-1}$
由于转移矩阵的特殊性,不难证明$a_i$恰好是线性递推公式里的系数。
假设我们已经将$M^n$表示为 $b_0 M^0\\ +\\ b_1 M^1\\ +\\ \\cdots\\ b_{K-1}M^{K-1}$这样的形式,不难得到$M^{n+1}$的表示法。只要将$M^n$乘个M之后得到的项中$M^K$拆成小于K次的线性组合就好了。 这样我们可以预处理出$M^0\\ M^1\\ \\cdots\\ M^{2K-2}$的表示法。
对于次数更高的, $M^{i+j}=M^i*M^j$ 可以看成是两个多项式的乘法。 利用快速幂 可以在O(K2logn)的时间求出$M^n$的表示法.
另外有一个优化常数的trick, 可以预处理出$M^1$ $M^2$ $M^4$ $M^8$.... $M^{2^r}$这些项, 对于$M^n$只要根据二进制位相应的乘上这些项就好了。 这样做比直接做快速幂快一倍(少了一半的多项式乘法操作)。
参考代码：
1 //ans=12747994
2 #include <cstdio>
3 #include <iostream>
4 #include <queue>
5 #include <algorithm>
6 #include <cstring>
7 #include <set>
8 using namespace std;
9
10 #define N 2000
11 typedef long long ll;
12
13 const int Mod=20092010;
14 int a[N],f[N<<1];
15 int k=2000;
16
17
18
19 //基本思想是把A^n 表示成A^0 A^1 A^2 ... A^(k-1)的线性组合
20 //A^(p+q)可看成两个多项式相乘,只要实现预处理出A^0 A^1 A^2 ... A^(2k-2)的多项式表示法
21 //A^k可以根据特征多项式的性质得到 ,A^(n+1)次可以从A^n次得到 根据这个来预处理
22 struct Poly
23 {
24 int b[N];
25 }P[N<<1];
26
27
28 Poly operator * (const Poly &A,const Poly &B)
29 {
30 Poly ans; memset(ans.b,0,sizeof(ans.b));
31 for (int i=0;i<=2*k-2;i++)
32 {
33 int res=0;
34 for (int j=max(0,i-k+1);j<k && j<=i;j++)
35 {
36 res+=1ll*A.b[j]*B.b[i-j]%Mod;
37 if (res>=Mod) res-=Mod;
38 }
39 if (i<k) {ans.b[i]=res; continue;}
40
41 //把次数大于等于k的搞成小于k
42 for (int j=0;j<k;j++)
43 {
44 ans.b[j]+=1ll*res*P[i].b[j]%Mod;
45 if (ans.b[j]>=Mod) ans.b[j]-=Mod;
46 }
47 }
48 return ans;
49 }
50
51 Poly Power_Poly(ll p)
52 {
53 if (p<=2*k-2) return P[p];
54
55 Poly ans=P[0],A=P[1];
56 for (;p;p>>=1)
57 {
58 if (p&1) ans=ans*A;
59 A=A*A;
60 }
61 return ans;
62 }
63
64 int main()
65 {
66 freopen("in.in","r",stdin);
67 freopen("out.out","w",stdout);
68
69 //f[n]=a[k-1]f[n-1]....a[0]f[n-k]
70 a[0]=a[1]=1; ll n; n=1e18;
71 for (int i=0;i<k;i++) P[i].b[i]=1;
72
73 //P[k]=a[0]P[0]+a[1]P[1]+....a[k-1]P[k-1]
74 for (int i=0;i<k;i++) P[k].b[i]=a[i];
75
76 //Calculate P[k+1]...P[2k-2]
77 //using P[n+1]=a[0]*b[k-1]+ (a[1]*b[k-1]+b[0]) + (a[2]*b[k-1]+b[1]) +...(a[k-1]*b[k-1]+b[k-2])
78 for (int j=k+1;j<=2*k-2;j++)
79 {
80 P[j].b[0]=1ll*a[0]*P[j-1].b[k-1]%Mod;
81 for (int i=1;i<k;i++)
82 P[j].b[i]=(1ll*a[i]*P[j-1].b[k-1]%Mod+P[j-1].b[i-1])%Mod;
83 }
84
85 Poly tmp=Power_Poly(n-k+1); int ans=0;
86
87 for (int i=0;i<k;i++) f[i]=1;
88 for (int i=k;i<=2*k-2;i++) f[i]=(f[i-1999]+f[i-2000])%Mod;
89
90 //A^n*X=b[0]*A^0*X+b[1]*A^1*X+...b[k-1]*A^(k-1)*X A^i*X= {f[k-1+i] f[k-2+i]... f[0+i]}
91 for (int i=0;i<k;i++)
92 {
93 ans+=1ll*tmp.b[i]*f[k-1+i]%Mod;
94 if (ans>=Mod) ans-=Mod;
95 }
96 printf("%d\\n",ans);
97 return 0;
98 }
2.设
，求
的值。其中
，
和
。
题目链接：http://www.51nod.com/onlineJudge/questionCode.html#!problemId=1229
首先这个题可以用扰动法 搞出一个关于k的递推公式,在O(k2)的时间解决。
具体可以参考这篇。 虽然不是同一个式子,但是方法是一样的,扰动法在《具体数学》上也有介绍。因此本文不再赘述。
给个AC代码供参考：
1 #include <cstdio> 2 #include <iostream> 3 #include <queue> 4 #include <algorithm> 5 #include <cstring> 6 #include <set> 7 using namespace std; 8 9 #define N 2010 10 typedef long long ll; 11 const int Mod=1000000007; 12 13 int k; 14 ll n,r; 15 int f[N],inv[N]; 16 int fac[N],fac_inv[N]; 17 18 int C(int x,int y) 19 { 20 if (y==0) return 1; 21 if (y>x) return 0; 22 23 int res=1ll*fac[x]*fac_inv[y]%Mod; 24 return 1ll*res*fac_inv[x-y]%Mod; 25 } 26 27 int Power(ll a,ll p) 28 { 29 int res=1; a%=Mod; 30 for (;p;p>>=1) 31 { 32 if (p&1) res=1ll*res*a%Mod; 33 a=a*a%Mod; 34 } 35 return res; 36 } 37 38 int Solve1() 39 { 40 f[0]=n; int t=n+1; 41 for (int i=1;i<=k;i++) 42 { 43 f[i]=t=1ll*t*(n+1)%Mod; 44 for (int j=0;j<i;j++) 45 { 46 f[i]+=Mod-1ll*C(i+1,j)*f[j]%Mod; 47 if (f[i]>=Mod) f[i]-=Mod; 48 } 49 f[i]--; if (f[i]<0) f[i]+=Mod; 50 f[i]=1ll*f[i]*inv[i+1]%Mod; 51 } 52 return f[k]; 53 } 54 55 56 int Solve2() 57 { 58 f[0]=Power(r,n+1)-r%Mod; 59 if (f[0]<0) f[0]+=Mod; 60 f[0]=1ll*f[0]*Power(r-1,Mod-2)%Mod; 61 62 for (int i=1;i<=k;i++) 63 { 64 f[i]=1ll*Power(n+1,i)*Power(r,n+1)%Mod; 65 f[i]-=r%Mod; if (f[i]<0) f[i]+=Mod; 66 67 int tmp=0; 68 for (int j=0;j<i;j++) 69 { 70 tmp+=1ll*C(i,j)*f[j]%Mod; 71 if (tmp>=Mod) tmp-=Mod; 72 } 73 f[i]-=1ll*(r%Mod)*tmp%Mod; 74 if (f[i]<0) f[i]+=Mod; 75 f[i]=1ll*f[i]*Power(r-1,Mod-2)%Mod; 76 //cout<<i<<" "<<f[i]<<endl; 77 } 78 return f[k]; 79 } 80 81 82 83 int main() 84 { 85 //freopen("in.in","r",stdin); 86 //freopen("out.out","w",stdout); 87 88 inv[1]=1; for (int i=2;i<N;i++) inv[i]=1ll*(Mod-Mod/i)*inv[Mod%i]%Mod; 89 fac[0]=1; for (int i=1;i<N;i++) fac[i]=1ll*fac[i-1]*i%Mod; 90 fac_inv[0]=1; for (int i=1;i<N;i++) fac_inv[i]=1ll*fac_inv[i-1]*inv[i]%Mod; 91 92 int T; scanf("%d",&T); 93 while (T--) 94 { 95 cin >> n >> k >> r; 96 if (r==1) printf("%d\\n",Solve1()); 97 else printf("%d\\n",Solve2()); 98 } 99 100 return 0; 101 }
本文要介绍的是利用优化后的矩阵乘法来解决本题(也许是我写的太丑,常数巨大,极限数据要跑6s,不能AC本题,但是可以拿来作为练习)
首先要知道怎么构造矩阵。
设$F(n,j)=n^j r^n$ 列向量 $X=[S_{n-1}\\ ,\\ F(n,k)\\ ,\\ F(n,k-1),\\ \\cdots\\ ,\\ F(n,0)]^T$
更加详细的题解可以参考这篇. 上面的推导过程的出处(矩阵实在不会用latex公式打,就只好copy了。)
我的方法和他略有不同, 我的列向量第一个元素是$S_{n-1}$ ,因此我的转移矩阵的第一行是1,1,0,0...0 其他都是一样的。
另外我猜测这题的这个式子和国王奇遇记那题一样,应该也是一个什么玩意乘上一个多项式,可以用多项式插值的办法来求(排行榜前面的代码跑的都很快,目测是O(K)的)。
比较懒,懒得去想了。。有兴趣的朋友可以去搞一搞？
我的代码（最后几个点TLE,用来练习 优化矩阵乘法）：
1 //http://www.51nod.com/onlineJudge/questionCode.html#!problemId=1229 2 #include <iostream> 3 #include <cstdio> 4 #include <cmath> 5 #include <cstring> 6 #include <algorithm> 7 #include <vector> 8 #include <map> 9 #include <cstdlib> 10 #include <set> 11 using namespace std; 12 13 #define X first 14 #define Y second 15 #define Mod 1000000007 16 #define N 2010 17 typedef long long ll; 18 typedef pair<int,int> pii; 19 20 int K; 21 ll n,r; 22 int fac[N],fac_inv[N],S[N]; 23 24 struct Poly 25 { 26 int cof[N]; 27 void print() 28 { 29 for (int i=0;i<K;i++) printf("%d ",cof[i]); 30 printf("\\n"); 31 } 32 }M[N<<1],R[62]; 33 34 int a[N],b[N]; 35 36 int C(int x,int y) 37 { 38 if (y==0) return 1; 39 40 int res=1ll*fac[x]*fac_inv[y]%Mod; 41 return 1ll*res*fac_inv[x-y]%Mod; 42 } 43 44 int Power(ll x,ll p) 45 { 46 int res=1; x%=Mod; 47 for (;p;p>>=1) 48 { 49 if (p&1) res=1ll*res*x%Mod; 50 x=x*x%Mod; 51 } 52 return res; 53 } 54 55 56 Poly operator * (const Poly &A,const Poly &B) 57 { 58 Poly ans; memset(ans.cof,0,sizeof(ans.cof)); 59 60 for (int i=0;i<2*K-1;i++) 61 { 62 int res=0; 63 for (int j=max(0,i-K+1);j<=i && j<K;j++) 64 { 65 res+=1ll*A.cof[j]*B.cof[i-j]%Mod; 66 if (res>=Mod) res-=Mod; 67 } 68 if (i<K) {ans.cof[i]=res; continue;} 69 70 for (int j=0;j<K;j++) 71 { 72 ans.cof[j]+=1ll*res*M[i].cof[j]%Mod; 73 if (ans.cof[j]>=Mod) ans.cof[j]-=Mod; 74 } 75 } 76 return ans; 77 } 78 79 Poly Poly_Power(ll p) 80 { 81 if (p<=2*K-2) return M[p]; 82 83 Poly res=M[0]; 84 //p&(1ll<<j) 千万别忘了1后面的ll ！！！ 85 for (int j=0;j<=60;j++) if (p&(1ll<<j)) res=res*R[j]; 86 return res; 87 } 88 89 90 void Init() 91 { 92 memset(a,0,sizeof(a)); 93 memset(b,0,sizeof(b)); 94 95 96 //求出特征多项式 M^(K+2)=a[0]*M^0 + a[1]*M^1 + ... a[K+1]*M^(K+1) 97 int op; 98 if (r==1) //特征多项式f(x)= (x-1)^(K+2) 99 { 100 for (int i=0;i<K+2;i++) 101 { 102 op=(K-i)&1? -1:1; 103 a[i]=op*C(K+2,i); 104 a[i]=-a[i]; 105 if (a[i]<0) a[i]+=Mod; 106 } 107 } 108 else //特征多项式f(x)= (x-1) * (x-r)^(K+1) 109 { 110 for (int i=0;i<=K+1;i++) 111 { 112 op=(K+1-i)&1? -1:1; 113 b[i]=1ll*op*C(K+1,i)*Power(r,K+1-i)%Mod; 114 if (b[i]<0) b[i]+=Mod; 115 } 116 for (int i=1;i<K+2;i++) a[i]=b[i-1]; 117 for (int i=0;i<=K+1;i++) 118 { 119 a[i]-=b[i]; 120 a[i]=-a[i]; 121 if (a[i]<0) a[i]+=Mod; 122 } 123 } 124 125 K+=2; //矩阵的规格 126 127 128 //预处理M^0...M^(2K-2) 129 130 memset(M,0,sizeof(M)); 131 132 for (int i=0;i<K;i++) M[i].cof[i]=1; 133 for (int i=0;i<K;i++) M[K].cof[i]=a[i]; 134 135 for (int i=K+1;i<=2*K-2;i++) 136 { 137 M[i].cof[0]=1ll*a[0]*M[i-1].cof[K-1]%Mod; 138 for (int j=1;j<K;j++) 139 { 如何更改theorem 环境James Munkres Topology: Theorem 20.4