常係數齊次線性遞推式第n項的快速計算初探 XJB學後的XBJ胡扯
要做啥？
求(f[n]=sum_{i=1}^{k} a[i]f[n-i] ;(nge k))，(a[1,k],f[0,k-1])已經給出。
我會矩陣快速冪！
時間複雜度(O(k^3log n))，其中(nle 10^9,kle32000)，emmm。
我會魔法！
設初始狀態矩陣為(S)，轉移矩陣為(A)。轉移長得像這個樣子（四階情形）
[
egin{bmatrix}
f[4]\
f[3]\
f[2]\
f[1]
end{bmatrix}
=egin{bmatrix}
a[1]&a[2]&a[3]&a[4]\
1&0&0&0\
0&1&0&0\
0&0&1&0\
end{bmatrix}
egin{bmatrix}
f[3]\
f[2]\
f[1]\
f[0]
end{bmatrix}
]
構造序列(c)滿足(A^n=sum_{i=0}^{k-1}c[i]A^i)，兩邊同時右乘(S)，
[
A^nS=(sum_{i=0}^{k-1}c[i]A^i)S=sum_{i=0}^{k-1}c[i]A^iS
]
我們的答案為((A^nS)[k-1,0])，注意(A)的特徵，可以發現，
[
(A^nS)[k-1,0]=sum_{i=0}^{k-1} c[i](A^iS)[k-1,0]=sum_{i=0}^{k-1}c[i]f[i]
]
所以只要構造出(c)，我們就能(O(k))地完成答案的計算。
設(k)次多項式(G(x))滿足(G(A)=O)，則有
[
R(x)=x^nmod G(x)\
R(A)=A^n
]
顯然(R(x))是個(k-1)次多項式，(c)即為(R(x))的係數。
即需要構造一個滿足要求的(G(x))。
我會更強力的魔法！
(k)階方陣(A)的特徵多項式 (p(lambda)=det(lambda I-A))，顯然(p(lambda))是一個(k)次多項式。
哈密爾頓–凱萊定理 (p(A)=O)。（可以舉例子簡單驗證一下）
於是(G(x))就這樣出來了，推式子可得
[
p(lambda)=lambda^k-sum_{i=1}^ka[i]lambda^{k-i}
]
然後回代就做完了。取模的那部分是快速冪套多項式取模，所以總的時間複雜度為(O(klog klog n))。
我不會魔法的實現
暫時留坑 題目傳送門
已填坑 on 26 Feb。參考程式碼如下
int n,k,p;
int a[N],f[N],c[N];
int g[N],invg[N];
void pomo(int*F,int*R) {
static int Q[N],T[N];
reverse(F,F+k+k-1);
copy(F,F+k,T);
nuthtr(T,p,1);
for(int i=0; i<p; ++i) Q[i]=1LL*T[i]*invg[i]%P;
nuthtr(Q,p,-1);
fill(Q+k-1,Q+p,0);
reverse(F,F+k+k-1);
reverse(Q,Q+k-1);
nuthtr(Q,p,1);
for(int i=0; i<p; ++i) Q[i]=1LL*g[i]*Q[i]%P;
nuthtr(Q,p,-1);
for(int i=0; i<k; ++i) R[i]=(F[i]-Q[i]+P)%P;
fill(R+k,R+p,0);
fill(Q,Q+p,0);
fill(T,T+p,0);
}
int main() {
n=read(),k=read();
for(int i=1; i<=k; ++i) a[i]=(read()%P+P)%P;
for(int i=0; i<k; ++i) f[i]=(read()%P+P)%P;
for(int i=1; i<=k; ++i) g[k-i]=(P-a[i]); g[k]=1;
for(p=1; p<=k; p<<=1); p<<=1;
reverse(g,g+k+1);
poin(g,invg,p);
fill(invg+k+1,invg+p,0);
reverse(g,g+k+1);
nuthtr(g,p,1);
nuthtr(invg,p,1);
static int s[N];
s[1]=1;
c[0]=1;
for(; n; n>>=1) {
if(n&1) {
nuthtr(c,p,1);
nuthtr(s,p,1);
for(int i=0; i<p; ++i) c[i]=1LL*c[i]*s[i]%P;
nuthtr(c,p,-1);
nuthtr(s,p,-1);
pomo(c,c);
}
nuthtr(s,p,1);
for(int i=0; i<p; ++i) s[i]=1LL*s[i]*s[i]%P;
nuthtr(s,p,-1);
pomo(s,s);
}
int ans=0;
for(int i=0; i<k; ++i) ans=(ans+1LL*c[i]*f[i])%P;
printf("%d
",ans);
return 0;
}