利用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 }

时间： 2024-10-20 13:22:43

利用Cayley-Hamilton theorem 优化矩阵线性递推的相关文章

hdu1588---Gauss Fibonacci(矩阵,线性递推)

根据题意:最后就是求f(b) + f(k + b) + f(2 * k + b) + -+ f((n-1) * k + b) 显然f(b) = A^b 其中A = 1 1 1 0 所以sum(n - 1) = A^b(E + A^ k + A ^(2 * k) + - + A ^((n - 1) * k) 设D = A^k sum(n-1) = A^b(E + D + D ^ 2 + - + D ^(n - 1)) 括号里的部分就可以二分递归求出来了 而单个矩阵就可以用矩阵快速幂求出来 /***

多校第九场：贪心+矩阵快速幂中间优化+线性递推&amp;线段树递推

HDU 4968 Improving the GPA 思路:贪心的搞吧!比赛的时候想了好久,然后才发现了点规律,然后乱搞1A. 因为贪心嘛!大的情况就是刚开始每个人的分数都是最大的最小值,即绩点4.0的最低分数85,然后最后一个数设为剩余的分数,然后如果小于60就从第一个分数补到这个分数来,然后最后一个分数还小于60,那就用第二个补--依次往下搞,那时我也不知道这样就搞出答案了,我还没证明这个对不对呢,哈哈. 小的情况:小的情况就是先假设每个人都是绩点最小的最大分数,即绩点2.0的最大分数69,

HDU 5863 cjj&#39;s string game ( 16年多校10 G 题、矩阵快速幂优化线性递推DP )

题目链接 题意 : 有种不同的字符,每种字符有无限个,要求用这k种字符构造两个长度为n的字符串a和b,使得a串和b串的最长公共部分长度恰为m,问方案数 分析 : 直觉是DP 不过当时看到 n 很大.但是 m 很小的时候 发现此题DP并不合适.于是想可能是某种组合数学的问题可以直接公式算 看到题解的我.恍然大悟.对于这种数据.可以考虑一下矩阵快速幂优化的DP 首先要想到线性递推的 DP 式子 最直观的想法就是 dp[i][j] = 到第 i 个位置为止.前面最长匹配长度为 j 的方案数 但是如果仔

矩阵乘法优化线性递推

矩阵乘法是线性代数中一块很重要的内容.矩阵乘法的定义很奇怪[1],但正是这种奇怪的性质,让矩阵乘法成为在除了线性代数和其衍生学科(还有诸如矩阵力学之类)外最广泛使用的关于矩阵变换的应用.(什么?FFT不属于矩阵变换吧...) 注: [1]: 矩阵乘法有另外的很多定义,如未说明,指的是中间不带符号的矩阵乘法,即一般矩阵乘积.另有 标量乘积(即所有数乘上一个固定的数),阿达马乘积等,没有那么诡异,但是在大多数问题的用途上也不大. 你不会矩阵乘法?没关系,下一篇会写到的 矩阵乘法的本质 矩阵乘法的本质

矩阵乘法递推的优化艺术

对于一个线性递推式,求它第项的值,通常的做法是先构造一个的矩阵,然后在时间内求出. 其实,由于这个矩阵的特殊性,可以将时间优化到.接下来我会以一个题目来讲解矩阵乘法递推的优化. 题目:http://www.51nod.com/onlineJudge/questionCode.html#!problemId=1229 题意:设,求的值.其中,和 . 前言:本题如果用普通的矩阵做法,很明显会TLE.那么我们要对这个特殊的矩阵进行时间上的优化. 分析:本题主要可用两种方法解决,分别是错位相减和矩阵乘法

HDU 4914 Linear recursive sequence(矩阵乘法递推的优化)

题解见X姐的论文 矩阵乘法递推的优化,只是mark一下.. HDU 4914 Linear recursive sequence(矩阵乘法递推的优化)

HDU - 6172:Array Challenge (BM线性递推)

题意:给出,三个函数,h,b,a,然后T次询问,每次给出n,求sqrt(an); 思路:不会推,但是感觉a应该是线性的,这个时候我们就可以用BM线性递推,自己求出前几项,然后放到模板里,就可以求了. 数据范围在1e15,1000组都可以秒过. 那么主要的问题就是得确保是线性的,而且得求出前几项. #include<bits/stdc++.h> using namespace std; #define rep(i,a,n) for (int i=a;i<n;i++) #define per

矩阵经典题目七：Warcraft III 守望者的烦恼（矩阵加速递推）

https://www.vijos.org/p/1067 很容易推出递推式f[n] = f[n-1]+f[n-2]+......+f[n-k]. 构造矩阵的方法:构造一个k*k的矩阵,其中右上角的(k-1)*(k-1)的矩阵是单位矩阵,第k行的每个数分别对应f[n-1],f[n-2],,f[n-k]的系数.然后构造一个k*1的矩阵,它的第i行代表f[i],是经过直接递推得到的.设ans[][]是第一个矩阵的n-k次幂乘上第二个矩阵,f[n]就是ans[k][1]. 注意:用__int64 #in

[HDOJ6172] Array Challenge（线性递推，黑科技）

题目链接:http://acm.hdu.edu.cn/showproblem.php?pid=6172 题意:给一堆东西,就是求个线性递推式,求第n项%1e9+7 杜教板真牛逼啊,线性递推式用某特征值相关的论文板,打表前几项丢进去就出结果了. 1 #include <bits/stdc++.h> 2 using namespace std; 3 4 typedef long long ll; 5 #define rep(i,a,n) for (ll i=a;i<n;i++) 6 #def