首页 > 代码库 > 利用Cayley-Hamilton theorem 优化矩阵线性递推
利用Cayley-Hamilton theorem 优化矩阵线性递推
平时有关线性递推的题,很多都可以利用矩阵乘法来解k决。 时间复杂度一般是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
另外我在高等代数的课本上找到了 证明(和维基百科里的第一种证明方法是一样的)
下面介绍几个 利用可以这个定理解决的题目：
1. project euler 258
显然可以用矩阵乘法来做。下面讲一下怎么利用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 { 140 M[i].cof[j]=1ll*a[j]*M[i-1].cof[K-1]%Mod+M[i-1].cof[j-1]; 141 if (M[i].cof[j]>=Mod) M[i].cof[j]-=Mod; 142 } 143 } 144 145 //预处理M^1 M^2 M^4 M^8 M^16... 146 R[0]=M[1]; 147 for (int i=1;i<=60;i++) R[i]=R[i-1]*R[i-1]; 148 149 S[0]=0; 150 for (int i=1;i<K;i++) 151 { 152 S[i]=1ll*Power(i,K-2)*Power(r,i)%Mod; 153 S[i]+=S[i-1]; 154 if (S[i]>=Mod) S[i]-=Mod; 155 } 156 } 157 158 159 int Solve() 160 { 161 Poly tmp=Poly_Power(n); 162 163 int ans=0,t; 164 165 for (int i=0;i<K;i++) 166 { 167 ans+=1ll*tmp.cof[i]*S[i]%Mod; 168 if (ans>=Mod) ans-=Mod; 169 } 170 return ans; 171 } 172 173 int main() 174 { 175 //freopen("in.in","r",stdin); 176 //freopen("out.out","w",stdout); 177 178 fac[0]=1; 179 for (int i=1;i<N;i++) fac[i]=1ll*fac[i-1]*i%Mod; 180 fac_inv[N-1]=Power(fac[N-1],Mod-2); 181 for (int i=N-2;i>=0;i--) fac_inv[i]=1ll*fac_inv[i+1]*(i+1)%Mod; 182 183 184 int T; scanf("%d",&T); 185 while (T--) 186 { 187 cin >> n >> K >> r; 188 Init(); 189 printf("%d\n",Solve()); 190 } 191 192 193 return 0; 194 }
3. hdu 3483
福利,和上面那题几乎完全是一样的,就是 范围小了很多。可以用来测试上面未通过的代码。 坑点就是 模数 最大是2e9, 最好开个long long, 否则int范围做加法就会爆。
参考代码：
1 #include <iostream>
2 #include <cstdio>
3 #include <cmath>
4 #include <cstring>
5 #include <algorithm>
6 #include <vector>
7 #include <map>
8 #include <cstdlib>
9 #include <set>
10 using namespace std;
11
12 #define X first
13 #define Y second
14 #define N 70
15 typedef long long ll;
16 typedef pair<int,int> pii;
17
18 int K,Mod;
19 ll n,r;
20 ll S[N];
21 ll cob[N][N];
22
23 struct Poly
24 {
25 ll cof[N];
26 void print()
27 {
28 for (int i=0;i<K;i++) printf("%d ",cof[i]);
29 printf("\n");
30 }
31 }M[N<<1],R[62];
32
33 ll a[N],b[N];
34
35 ll C(int x,int y)
36 {
37 return cob[x][y];
38 }
39
40 ll Power(ll x,ll p)
41 {
42 ll res=1; x%=Mod;
43 for (;p;p>>=1)
44 {
45 if (p&1) res=1ll*res*x%Mod;
46 x=x*x%Mod;
47 }
48 return res;
49 }
50
51
52 Poly operator * (const Poly &A,const Poly &B)
53 {
54 Poly ans; memset(ans.cof,0,sizeof(ans.cof));
55
56 for (int i=0;i<2*K-1;i++)
57 {
58 ll res=0;
59 for (int j=max(0,i-K+1);j<=i && j<K;j++)
60 {
61 res+=1ll*A.cof[j]*B.cof[i-j]%Mod;
62 if (res>=Mod) res-=Mod;
63 }
64 if (i<K) {ans.cof[i]=res; continue;}
65
66 for (int j=0;j<K;j++)
67 {
68 ans.cof[j]+=1ll*res*M[i].cof[j]%Mod;
69 if (ans.cof[j]>=Mod) ans.cof[j]-=Mod;
70 }
71 }
72 return ans;
73 }
74
75 Poly Poly_Power(ll p)
76 {
77 if (p<=2*K-2) return M[p];
78
79 Poly res=M[0];
80 //p&(1ll<<j) 千万别忘了1后面的ll ！！！
81 for (int j=0;j<=60;j++) if (p&(1ll<<j)) res=res*R[j];
82 return res;
83 }
84
85
86 void Init()
87 {
88 memset(a,0,sizeof(a));
89 memset(b,0,sizeof(b));
90
91 cob[0][0]=1;
92 for (int i=1;i<N;i++)
93 {
94 cob[i][0]=1;
95 for (int j=1;j<N;j++)
96 cob[i][j]=(cob[i-1][j-1]+cob[i-1][j])%Mod;
97 }
98
99 //求出特征多项式 M^(K+2)=a[0]*M^0 + a[1]*M^1 + ... a[K+1]*M^(K+1)
100 int op;
101 if (r==1) //特征多项式f(x)= (x-1)^(K+2)
102 {
103 for (int i=0;i<K+2;i++)
104 {
105 op=(K-i)&1? -1:1;
106 a[i]=op*C(K+2,i);
107 a[i]=-a[i];
108 if (a[i]<0) a[i]+=Mod;
109 }
110 }
111 else //特征多项式f(x)= (x-1) * (x-r)^(K+1)
112 {
113 for (int i=0;i<=K+1;i++)
114 {
115 op=(K+1-i)&1? -1:1;
116 b[i]=1ll*op*C(K+1,i)*Power(r,K+1-i)%Mod;
117 if (b[i]<0) b[i]+=Mod;
118 }
119 for (int i=1;i<K+2;i++) a[i]=b[i-1];
120 for (int i=0;i<=K+1;i++)
121 {
122 a[i]-=b[i];
123 a[i]=-a[i];
124 if (a[i]<0) a[i]+=Mod;
125 }
126 }
127
128 K+=2; //矩阵的规格
129
130
131 //预处理M^0...M^(2K-2)
132
133 memset(M,0,sizeof(M));
134
135 for (int i=0;i<K;i++) M[i].cof[i]=1;
136 for (int i=0;i<K;i++) M[K].cof[i]=a[i];
137
138 for (int i=K+1;i<=2*K-2;i++)
139 {
140 M[i].cof[0]=1ll*a[0]*M[i-1].cof[K-1]%Mod;
141 for (int j=1;j<K;j++)
142 {
143 M[i].cof[j]=1ll*a[j]*M[i-1].cof[K-1]%Mod+M[i-1].cof[j-1];
144 if (M[i].cof[j]>=Mod) M[i].cof[j]-=Mod;
145 }
146 }
147
148 //预处理M^1 M^2 M^4 M^8 M^16...
149 R[0]=M[1];
150 for (int i=1;i<=60;i++) R[i]=R[i-1]*R[i-1];
151
152 S[0]=0;
153 for (int i=1;i<K;i++)
154 {
155 S[i]=1ll*Power(i,K-2)*Power(r,i)%Mod;
156 S[i]+=S[i-1];
157 if (S[i]>=Mod) S[i]-=Mod;
158 }
159 }
160
161
162 ll Solve()
163 {
164 Poly tmp=Poly_Power(n);
165
166 ll ans=0;
167
168 for (int i=0;i<K;i++)
169 {
170 ans+=1ll*tmp.cof[i]*S[i]%Mod;
171 if (ans>=Mod) ans-=Mod;
172 }
173 return ans;
174 }
175
176 int main()
177 {
178 //freopen("in.in","r",stdin);
179 //freopen("out.out","w",stdout);
180
181 int T;
182 while (true)
183 {
184 cin >> n >> K >> Mod; r=K;
185 if (n<0) break;
186 Init();
187 printf("%I64d\n",Solve());
188 }
189
190
191 return 0;
192 }
4.codechef Dec Challenge POWSUMS https://www.codechef.com/problems/POWSUMS
官方题解： https://discuss.codechef.com/questions/86250/powsums-editorial
参考代码：
1 //https://www.codechef.com/problems/POWSUMS
2 //https://discuss.codechef.com/questions/86250/powsums-editorial
3 //https://discuss.codechef.com/questions/49614/linear-recurrence-using-cayley-hamilton-theorem
4
5
6 #include <iostream>
7 #include <cstdio>
8 #include <cmath>
9 #include <cstring>
10 #include <algorithm>
11 #include <vector>
12 #include <map>
13 #include <cstdlib>
14 #include <set>
15 using namespace std;
16
17 #define X first
18 #define Y second
19 #define Mod 1000000007
20 #define N 310
21 typedef long long ll;
22 typedef pair<int,int> pii;
23
24 inline int Mul(int x,int y){return 1ll*x*y%Mod;}
25 inline int Add(int x,int y){return ((x+y)%Mod+Mod)%Mod;}
26
27 int n,Q;
28 int a[N],e[N],inv[N],f[N<<1];
29
30 struct Poly
31 {
32 int cof[N];
33 void print()
34 {
35 for (int i=0;i<n;i++) printf("%d ",cof[i]);
36 printf("\n");
37 }
38 }M[N<<1],tt[63];
39
40 Poly operator * (const Poly &A,const Poly &B)
41 {
42 Poly ans; memset(ans.cof,0,sizeof(ans.cof));
43
44 for (int i=0;i<=2*n-2;i++)
45 {
46 int res=0;
47 for (int j=max(0,i-n+1);j<n && j<=i;j++)
48 {
49 res+=1ll*A.cof[j]*B.cof[i-j]%Mod;
50 if (res>=Mod) res-=Mod;
51 }
52 if (i<n) {ans.cof[i]=res;continue;}
53 for (int j=0;j<n;j++)
54 {
55 ans.cof[j]+=1ll*res*M[i].cof[j]%Mod;
56 if (ans.cof[j]>=Mod) ans.cof[j]-=Mod;
57 }
58 }
59 return ans;
60 }
61
62 void Init()
63 {
64 scanf("%d%d",&n,&Q);
65 for (int i=1;i<=n;i++) scanf("%d",&f[i]);
66
67 e[0]=1;
68 for (int i=1;i<=n;i++)
69 {
70 int op=1; e[i]=0;
71 for (int j=1;j<=i;j++,op=-op)
72 {
73 e[i]+=1ll*op*e[i-j]*f[j]%Mod;
74 if (e[i]<0) e[i]+=Mod;
75 if (e[i]>=Mod) e[i]-=Mod;
76 }
77 e[i]=1ll*e[i]*inv[i]%Mod;
78 //cout<<i<<" "<<e[i]<<endl;
79 }
80
81 for (int i=n+1;i<=2*n-1;i++)
82 {
83 f[i]=0; int op=1;
84 for (int j=1;j<=n;j++,op=-op)
85 {
86 f[i]+=1ll*op*e[j]*f[i-j]%Mod;
87 if (f[i]>=Mod) f[i]-=Mod;
88 if (f[i]<0) f[i]+=Mod;
89 }
90 }
91
92 for (int i=0;i<n;i++)
93 for (int j=0;j<n;j++)
94 M[i].cof[j]= (i==j);
95
96 //M^n= sum (a[i]*A^i)
97 for (int i=n-1,op=1;i>=0;i--,op=-op)
98 {
99 int tmp=op*e[n-i];
100 if (tmp<0) tmp+=Mod;
101 M[n].cof[i]=a[i]=tmp;
102 }
103
104 //Calc linear combination form of M^(n+1)...M^(2n-2)
105 for (int i=n+1;i<=2*n-2;i++)
106 {
107 M[i].cof[0]=1ll*a[0]*M[i-1].cof[n-1]%Mod;
108 for (int j=1;j<n;j++)
109 {
110 M[i].cof[j]=1ll*a[j]*M[i-1].cof[n-1]%Mod+M[i-1].cof[j-1];
111 if (M[i].cof[j]>=Mod) M[i].cof[j]-=Mod;
112 }
113
114 }
115
116 //预处理出M^2 M^4 M^8... 随机数据可以加速很多
117 tt[0]=M[1];
118 for (int i=1;i<=60;i++) tt[i]=tt[i-1]*tt[i-1];
119 }
120
121 Poly Poly_Power(ll p)
122 {
123 if (p<=2*n-2) return M[p];
124 Poly res=M[0];
125 for (int j=60;j>=0;j--) if (p&(1ll<<j)) res=res*tt[j];
126 return res;
127 }
128
129 int Solve(ll x)
130 {
131 if (x<=2*n-2) return f[x];
132
133 Poly tmp=Poly_Power(x-n);
134
135 int res=0;
136 for (int i=0;i<n;i++)
137 {
138 res+=1ll*tmp.cof[i]*f[n+i]%Mod;
139 if (res>=Mod) res-=Mod;
140 }
141
142 if (res<0) res=res%Mod+Mod;
143
144 return res;
145 }
146
147 int main()
148 {
149 //freopen("in.in","r",stdin);
150 //freopen("out.out","w",stdout);
151
152 //calc_inv
153 inv[1]=1;
154 for (int i=2;i<N;i++) inv[i]=1ll*(Mod-Mod/i)*inv[Mod%i]%Mod;
155 int T; ll x;
156
157 scanf("%d",&T);
158 while (T--)
159 {
160 Init();
161 while (Q--)
162 {
163 scanf("%lld",&x);
164 printf("%d ",Solve(x));
165 }
166 printf("\n");
167 }
168
169 return 0;
170 }
利用Cayley-Hamilton theorem 优化矩阵线性递推
声明：以上内容来自用户投稿及互联网公开渠道收集整理发布，本网站不拥有所有权，未作人工编辑处理，也不承担相关法律责任，若内容有误或涉及侵权可进行投诉： 投诉/举报 工作人员会在5个工作日内联系你，一经查实，本站将立刻删除涉嫌侵权内容。