@description@
给定递推关系式：
并给定 (A_1, A_2, dots , A_k) 的值，求 (A_n) 的值模 104857601。
input：
第一行给出两个整数 n，k。
第二行包含 k 个整数 A1，A2，...，Ak。
第三行包含 k 个整数 C1，C2，...，Ck。
output：
输出 An 的值。
constraints：
1 ≤ n ≤ 10^18
0 ≤ Ai, Ci < 104857601
1 ≤ k ≤ 30000
sample input：
3 5
1 2 3
4 5 6
sample output：
139
sample explanation：
(A_1 = 1,A2 = 2,A3 = 3)
(A4 = (3 × 4 + 2 × 5 + 1 × 6) mod 104857601 = 28)
(A5 = (28 × 4 + 3 × 5 + 2 × 6) mod 104857601 = 139)
@solution@
这道题，如果你学过常系数齐次线性递推，就是一道模板题。
所以以下内容我就来 BB 什么是常系数齐次线性递推。
@part - 1@
解决这道题的经典算法就是矩阵乘法（也可以分治 FFT 或者多项式求逆）。
设出系数矩阵
并设 (B_i = egin{bmatrix}A_{i+k-1}\A_{i+k-2}\vdots\A_{i+1}\A_{i}\end{bmatrix})
则：
即 (MB_{i-k}=B_{i-k+1})，或者等价地 (MB_{i}=B_{i+1})
于是就有 (M^{n-1}B_1=B_n)，而 (B_1) 是已知的，做一个矩阵快速幂求 (M^{n-1}) 即可。
然而这样的复杂度为 (O(k^3log n))，对于这道题并不能通过。
我们接下来就来优化这个算法。
@part - 2@
【注：下面多项式中 x 并不是代表一个实数，而是一个什么都可以代表的东西。我一开始并没有理解，所以懵逼了好久 QAQ……】
定义多项式：
通过一系列证明（证明我在下面part 3里面弄），发现（右边那个 “0” 不是个数字，而是表示零矩阵）：
我们用多项式除法，用 (x^n) 除以 (f(x)) 得到商式 (g(x)) 与余式 (r(x))。即：
代入矩阵 M 得到：
又因为 (f(M) = 0)，所以得到 (M^n=r(M))。
其中多项式 (r(x)) 的最高次数严格小于 k。
然后我们就把算法优化到 (O(k^4)) 的了……？
我们再来看看我们要求解的是什么（下面那个 (r(M)) 是 (M^{n-1}) 的余式）：
我们只需要 (B_n) 中第最后一行的 (A_n)，所以有：
时间复杂度在于多项式快速幂和取模 (O(klog klog n)) 。
暴力写取模 (O(k^2log n)) ，以牺牲时间为代价降低代码复杂度，但是对于这道题而言是过不了的。
一个优化矩阵的算法最后代码全程不见矩阵……
@part - 3@
本节属于理性愉悦。
先给出几个概念吧：
特征值：若常数 (lambda) 对于矩阵 (M)，存在向量 (vec x) 使得：
则称 (lambda) 为矩阵 (M) 的特征值。
变形上式得 ((lambda I - M)x = 0)，该方程有解的充要条件为 (det(lambda I - M) = 0)。
特征多项式：(det(lambda I - M)) 其实是一个多项式，我们称这个多项式为矩阵 M 的特征多项式，记为 (p(x))。
Cayley - Hamilton 定理：
浅显，但并不易懂。该定理的证明超出了我的理解范围。所以这里并没有证明。
大家的数学知识如果比我不知道高到哪里去了，可以看这一篇博客。
update in 2020/07/17：大概理解了一下百度百科的第一个证明，发现其实就是利用伴随矩阵的性质列等式，对比左右系数然后暴算。。。
上面提到的那位博主更新了一篇新博客，里面有这样一段话：
为什么会想到这么证明呢？是因为我们要把特征多项式和矩阵乘积连接起来，然后用矩阵乘积的方法来证明。 而伴随矩阵恰好可以把行列式计算转化为矩阵乘积。观察以上证明过程 (B) 的作用主要是连接 ((λI−A)) 和 (∣λI−A∣)。
由于矩阵的本质是线性变换，所以也可以从线性变换的角度理解。不过我线代太差了所以不太懂。。。
对于矩阵 (M)：
如果要求解它的特征多项式，则：
以第一行展开求行列式，得到：
其中 (M_{ij}) 是 (M) 的代数余子式。
可以发现去掉第一行和第 i 列，剩下的矩阵是一个下三角矩阵。直接对角线相乘就可以求出行列式了。
所以：
结合上面的 Cayley - Hamilton 定理，我们就得到了我们想要的东西。
@accepted code@
#include<cstdio>
#include<algorithm>
using namespace std;
typedef long long ll;
const int G = 3;
const int MAXN = 30000*5;
const int MOD = 104857601;
int pow_mod(int b, int p) {
int ret = 1;
while( p ) {
if( p & 1 ) ret = 1LL*ret*b%MOD;
b = 1LL*b*b%MOD;
p >>= 1;
}
return ret;
}
struct Polynomial{
void poly_copy(int *A, int *B, int n) {
for(int i=0;i<n;i++)
A[i] = B[i];
}
void poly_clear(int *A, int l, int r) {
for(int i=l;i<r;i++)
A[i] = 0;
}
void poly_revcopy(int *A, int *B, int n) {
for(int i=0;i<n;i++)
A[i] = B[n-i-1];
}
void ntt(int *A, int n, int type) {
for(int i=0,j=0;i<n;i++) {
if( i < j ) swap(A[i], A[j]);
for(int l=(n>>1);(j^=l)<l;l>>=1);
}
for(int s=2;s<=n;s<<=1) {
int t = (s>>1);
int u = (type == -1) ? pow_mod(G, (MOD-1) - (MOD-1)/s) : pow_mod(G, (MOD-1)/s);
for(int i=0;i<n;i+=s) {
for(int j=0,p=1;j<t;j++,p=1LL*p*u%MOD) {
int x = A[i+j], y = 1LL*p*A[i+j+t]%MOD;
A[i+j] = (x + y)%MOD, A[i+j+t] = (x + MOD - y)%MOD;
}
}
}
if( type == -1 ) {
int inv = pow_mod(n, MOD-2);
for(int i=0;i<n;i++)
A[i] = 1LL*A[i]*inv%MOD;
}
}
int tmp1[MAXN + 5], tmp2[MAXN + 5];
void poly_mul(int *A, int *B, int *C, int n, int m) {
int len; for(len = 1;len < n+m-1;len <<= 1);
poly_copy(tmp1, A, n); poly_clear(tmp1, n, len);
poly_copy(tmp2, B, m); poly_clear(tmp2, m, len);
ntt(tmp1, len, 1); ntt(tmp2, len, 1);
for(int i=0;i<len;i++)
C[i] = 1LL*tmp1[i]*tmp2[i]%MOD;
ntt(C, len, -1);
}
int tmp3[MAXN + 5];
void poly_inv(int *A, int *B, int n) {
if( n == 1 ) {
B[0] = pow_mod(A[0], MOD-2);
return ;
}
int len; for(len = 1;len < (n<<1);len <<= 1);
poly_inv(A, B, (n + 1) >> 1);
poly_copy(tmp3, A, n); poly_clear(tmp3, n, len);
ntt(tmp3, len, 1); ntt(B, len, 1);
for(int i=0;i<len;i++) B[i] = 1LL*B[i]*(2 + MOD - 1LL*tmp3[i]*B[i]%MOD)%MOD;
ntt(B, len, -1); poly_clear(B, n, len);
}
int tmp4[MAXN + 5], tmp5[MAXN + 5], tmp6[MAXN + 5];
void poly_mod(int *A, int *B, int *R, int n, int m) {
poly_revcopy(tmp4, B, m); poly_clear(tmp5, 0, 2*(n-m+1)); poly_inv(tmp4, tmp5, n-m+1);
poly_revcopy(tmp4, A, n); poly_mul(tmp4, tmp5, tmp6, n-m+1, n-m+1);
poly_revcopy(tmp4, tmp6, n-m+1); poly_copy(tmp5, B, m);
poly_mul(tmp4, tmp5, tmp6, n-m+1, m);
for(int i=0;i<m-1;i++)
R[i] = (A[i] + MOD - tmp6[i])%MOD;
}
int tmp7[MAXN + 5], tmp8[MAXN + 5];
void poly_pow(int *A, int *B, int *mod, ll p, int n) {
poly_copy(tmp7, A, n); B[0] = 1;
while( p ) {
if( p & 1 ) {
poly_mul(B, tmp7, tmp8, n, n), poly_copy(B, tmp8, 2*n-1);
poly_mod(B, mod, tmp8, 2*n-1, n), poly_copy(B, tmp8, n-1), poly_clear(B, n-1, 2*n-1);
}
poly_mul(tmp7, tmp7, tmp8, n, n), poly_copy(tmp7, tmp8, 2*n-1);
poly_mod(tmp7, mod, tmp8, 2*n-1, n), poly_copy(tmp7, tmp8, n-1), poly_clear(tmp7, n-1, 2*n-1);
p >>= 1;
}
}
}poly;
int A[MAXN + 5], C[MAXN + 5];
int f[MAXN + 5], p[MAXN + 5], r[MAXN + 5];
int main() {
ll N; int K;
scanf("%d%lld", &K, &N);
for(int i=1;i<=K;i++)
scanf("%d", &A[i]);
for(int i=1;i<=K;i++)
scanf("%d", &C[i]);
for(int i=0;i<K;i++)
f[i] = (MOD - C[K-i]);
f[K] = p[1] = 1;
poly.poly_pow(p, r, f, N-1, K+1);
int ans = 0;
for(int i=0;i<K;i++)
ans = (ans + 1LL*r[i]*A[i+1]%MOD)%MOD;
printf("%d
", ans);
}
@details@
注意因为要多次使用某些数组，所以用完就清是个好习惯（然后就被卡常）。
余式的次数是严格小于除式的次数，不能取等。