已知數列 $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 }

【U21729】遞推數列