已知数列 $f$ 满足
$$ f(0) = 0, f(1) = 1, f(n) = f(n-1) + 3f(n-2) $$
有 $t = {10} ^ 5$ 组询问，每次给定非负整数 $n \le {10} ^ {18}$ ，请你求出
$$
\sum_{i = 0} ^ {n-1} i f(i) \mod 10009
$$
做法 1 ：
$$\begin{pmatrix} & & 1 & & \\ & & & 1 & \\ 3 & & 1 & & \\ 6 & 3 & 1 & 1 & \\ & & & 1 & 1 \end{pmatrix} \begin{pmatrix} f_{n-2} \\ (n-2)f_{n-2} \\ f_{n-1} \\ (n-1)f_{n-1} \\ \sum_{i \le {n - 2} i f_i} \end{pmatrix} = \begin{pmatrix} f_{n-1} \\ (n-1)f_{n-1} \\ f_{n} \\ nf_{n} \\ \sum_{i \le {n-1} i f_i} \end{pmatrix}$$
复杂度为 $O(t \log n)$ ，带有 $125$ 倍的常数。
如果不使用 Caley-Hamition 定理，那么能不能过，就依赖于你的常数了。
矩阵乘法可以处理系数。
如果要处理 $\sum i ^ k f(i)$ ，考虑维护所有 $j \le k, \sum i ^ j f(i)$ ，根据二项式定理即可转移。
如果要处理 $\sum i ^ {\underline{k}} f(i)$ ，考虑描述成组合数的形式，根据二项式的加法公式进行转移。
1 #include <bits/stdc++.h> 2 using namespace std; 3 #define F(i, s, t) for (int i = (s), _ = (t); i <= _; i ++) 4 #define Fo(i, s, t) for (int i = (s), _ = (t); i < _; i ++) 5 #define LL long long 6 const int P = 10009; 7 struct Mat { 8 int h, w, a[5][5]; 9 Mat(int _1 = 0, int _2 = 0) : h(_1), w(_2) { memset(a, 0, sizeof a); } 10 int *operator [] (int x) { return a[x]; } 11 Mat operator * (Mat &B) { 12 Mat C(h, B.w); 13 Fo(k, 0, w) 14 Fo(i, 0, C.h) if (a[i][k]) 15 Fo(j, 0, C.w) 16 C.a[i][j] = (C.a[i][j] + (LL)a[i][k] * B.a[k][j]) % P; 17 return C; 18 } 19 }E, B, C, f; 20 namespace Input { 21 const int S = 10000000; 22 char s[S], *h = s+S, *t = h; 23 char nchar() { 24 if (h == t) fread(s, 1, S, stdin), h = s; 25 return *h ++; 26 } 27 LL rd() { 28 LL x = 0, f = 1; char c; 29 do { c = nchar(); if (c == ‘-‘) f = -1; } while (!isdigit(c)); 30 do { x = x*10+c-‘0‘, c = nchar(); } while (isdigit(c)); 31 return x * f; 32 } 33 } 34 using Input::rd; 35 Mat Pow(Mat A, LL n) { 36 if (! n) return E; 37 Mat t = Pow(A, n >> 1); 38 t = t * t; 39 if (n & 1) t = t * B; 40 return t; 41 } 42 int main() { 43 #ifndef ONLINE_JUDGE 44 freopen("729.in", "r", stdin); 45 #endif 46 E = Mat(5, 5), E[0][0] = E[1][1] = E[2][2] = E[3][3] = E[4][4] = 1; 47 B = Mat(5, 5), B[0][2] = B[1][3] = B[2][2] = B[3][2] = B[3][3] = B[4][3] = B[4][4] = 1, B[2][0] = B[3][1] = 3, B[3][0] = 6; 48 for (int cas = rd(); cas --; ) { 49 LL n = rd(); 50 if (! n) { puts("0"); continue; } 51 C = Pow(B, n-1); 52 f = Mat(5, 1), f[0][0] = f[1][0] = f[4][0] = 0, f[2][0] = f[3][0] = 1; 53 f = C * f; 54 printf("%d\n", f[4][0]); 55 } 56 return 0; 57 }
做法 2 ：
利用特征方程，解出 $f_n = \frac{1}{\sqrt {13}} (\frac{1 + \sqrt{13}}{2}) ^ n - \frac{1}{\sqrt{13}} (\frac{1 - \sqrt{13}}{2}) ^ n$ ，恰好 $13$ 在模 $10009$ 下有平方根，所以现在的问题变为求：
$$F(n) = \sum_{i = 0} ^ {n - 1} i q ^ i$$
一种做法是利用国王奇遇记的方法，一定存在一次多项式 $g(n)$ ，使得 $F(n) = q ^ n g(n) - g(0)$ ，求出 $g$ 的前两项然后插值，或者求出系数然后直接代入。
另外一种做法是求出 $i q ^ i$ 的前缀和（也就是逆差分函数），求的时候需要使用分部求和的技巧，
$$\Delta (uv) = Eu Ev - uv = Eu \Delta v + \Delta u v$$
$$uv = \Sigma Eu \Delta v + \Sigma \Delta uv$$
$$\Sigma (\Delta u) v = uv - \Sigma (Eu) (\Delta v)$$
$$\Sigma i q ^ i = \frac{q ^ i i}{q-1} - \Sigma \frac{q ^ {i + 1}}{q - 1}$$
$$\sum_{i = 0} ^ {n - 1} i q ^ i = \frac{q ^ n n}{q - 1} - \frac{q ^ {n + 1}}{(q - 1) ^ 2} + \frac{q}{(q - 1) ^ 2}$$
代入 $q = 2$ 进行检验是没错的，所以应该能够信任这个结果。
1 // f(n) = 4683 * 413 ^ n + -4683 * -412 ^ n 2 #include <bits/stdc++.h> 3 using namespace std; 4 #define LL long long 5 const int P = 10009; 6 namespace Input { 7 const int S = 10000000; 8 char s[S], *h = s+S, *t = h; 9 char nchar() { 10 if (h == t) fread(s, 1, S, stdin), h = s; 11 return *h ++; 12 } 13 LL rd() { 14 LL x = 0, f = 1; char c; 15 do { c = nchar(); if (c == ‘-‘) f = -1; } while (!isdigit(c)); 16 do { x = x*10+c-‘0‘, c = nchar(); } while (isdigit(c)); 17 return x * f; 18 } 19 } 20 using Input::rd; 21 int Pow(int x, LL t) { 22 int m = 1; 23 for (; t; t >>= 1, x = (LL)x * x % P) 24 if (t & 1) 25 m = (LL)m * x % P; 26 return m; 27 } 28 int Inv(int x) { return Pow(x, P-2); } 29 int C(LL n, int q) { 30 int I = Inv(q-1); 31 int u = Pow(q, n); 32 return ((LL)u * (n % P) * I - (LL)u * q * I * I + (LL)q * I * I) % P; 33 } 34 int T(LL n) { 35 int s = C(n, 413); 36 s = 4683 * (s - C(n, -412)) % P; 37 return s += (s < 0 ? P : 0); 38 } 39 int main() { 40 #ifndef ONLINE_JUDGE 41 freopen("729.in", "r", stdin); 42 #endif 43 for (int cas = rd(); cas --; ) { 44 LL n = rd(); 45 printf("%d\n", T(n)); 46 } 47 return 0; 48 }
