给 $n \times n(n \le 200000)$ 的方格，第一列、第一行 $f[1, i], f[i, 1] (\le 10 ^ 6)$ 和 $a, b, c(\le 10 ^ 6)$ ，转移方式为 $f[i, j] = a f[i, j-1] + b f[i-1, j] + c$ ，求 $f[n, n] \mod 1000003$ 。
考虑每个位置的贡献。
第一行、第一列的贡献直接 $O(n)$ 算，注意第一步的方向是确定的。
对于其余位置 $(i, j)$ ，贡献为 $\begin{aligned} c a ^ {n - j} b ^ {n - i} \binom{n - i + n - j}{n - i} \end{aligned}$ 。
$\begin{aligned} contr & = \sum_{i = 2} ^ {n} \sum_{j = 2} ^ n a ^ {n - j} b ^ {n - i} \binom{n - i + n - j}{n - i} \\ & = \sum_{i = 0} ^ {n - 2} \sum_{j = 0} ^ {n - 2} \binom{i + j}{i} a ^ i b ^ j \end{aligned}$
直观的想法是直接求 $\left\{ a_x \right\}$ 和 $\left\{ b_x \right\}$ 的二项卷积，系数相加就可以得到贡献。
实现的时候不能直接写 FFT ，因为会爆精度，需要使用「任意模数卷积」的技巧。
实现：
1 #include <bits/stdc++.h>
2 using namespace std;
3 #define F(i, a, b) for (register int i = (a), _b = (b); i <= _b; i++)
4 #define For(i, a, b) for (register int i = (a), _b = (b); i < _b; i++)
5 #define LL long long
6 #define db double
7 inline int rd(void) {
8 int f = 1; char c = getchar(); for (; !isdigit(c); c = getchar()) if (c == ‘-‘) f = -1;
9 int x = 0; for (; isdigit(c); c = getchar()) x = x*10+c-‘0‘; return x*f;
10 }
11
12 const int S = 600000;
13 const int P = 1e6+3;
14 const int p = sqrt(P);
15 const db PI = acos(-1);
16
17 int n, ans;
18 int fac[S], inv[S], ifac[S];
19 int m, a[S], b[S], I, c[S];
20 struct comp {
21 db x, y;
22 inline comp(db _x = 0, db _y = 0): x(_x), y(_y) {}
23 friend inline comp operator + (comp a, comp b) { return comp(a.x + b.x, a.y + b.y); }
24 friend inline comp operator - (comp a, comp b) { return comp(a.x - b.x, a.y - b.y); }
25 friend inline comp operator * (comp a, comp b) { return comp(a.x * b.x - a.y * b.y, a.x * b.y + a.y * b.x); }
26 }x[S], y[S], x2[S], y2[S], r[S];
27 int bit, len, rev[S];
28
29 inline int C(int n, int m) {
30 if (m < 0 || n < m) return 0;
31 return 1LL * fac[n] * ifac[m] % P * ifac[n-m] % P;
32 }
33
34 void FFT(comp *a, int sig) {
35 For(i, 0, len) if (i < rev[i]) swap(a[i], a[rev[i]]);
36 for (int i = 2; i <= len; i <<= 1) {
37 comp wn;
38 if (sig == 0)
39 wn = comp(cos(2*PI / +i), sin(2*PI / +i));
40 else wn = comp(cos(2*PI / -i), sin(2*PI / -i));
41 for (int j = 0; j < len; j += i) {
42 comp w(1, 0);
43 for (int k = 0; k < i / 2; k++) {
44 comp x = a[j+k], y = w * a[j+k+i/2];
45 a[j+k] = x+y, a[j+k+i/2] = x-y;
46 w = w * wn;
47 }
48 }
49 }
50 if (sig == 1) {
51 For(i, 0, len)
52 a[i] = comp(a[i].x / len, 0);
53 }
54 }
55 void contri(int t) {
56 For(i, 0, len) {
57 int w = (LL)(r[i].x + 0.5) % P;
58 c[i] = (c[i] + 1LL * w * t) % P;
59 }
60 }
61
62 int main(void) {
63 #ifndef ONLINE_JUDGE
64 freopen("4451.in", "r", stdin);
65 #endif
66
67 fac[0] = 1;
68 For(i, 1, S) fac[i] = 1LL * fac[i-1] * i % P;
69 inv[1] = 1;
70 For(i, 2, S) inv[i] = 1LL * (P - P / i) * inv[P % i] % P;
71 ifac[0] = 1;
72 For(i, 1, S) ifac[i] = 1LL * ifac[i-1] * inv[i] % P;
73
74 n = rd(), m = n-2, a[1] = rd(), b[1] = rd(), I = rd();
75 a[0] = 1; F(i, 2, n) a[i] = 1LL * a[i-1] * a[1] % P;
76 b[0] = 1; F(i, 2, n) b[i] = 1LL * b[i-1] * b[1] % P;
77
78 F(i, 1, n) {
79 int w = rd();
80 if (i > 1) ans = (ans + 1LL * w * a[n-1] % P * b[n-i] % P * C(n+n-i-2, n-2)) % P;
81 }
82 F(i, 1, n) {
83 int w = rd();
84 if (i > 1) ans = (ans + 1LL * w * a[n-i] % P * b[n-1] % P * C(n+n-i-2, n-2)) % P;
85 }
86
87 for (bit = 0, len = 1; len <= m+m; bit++, len <<= 1);
88 For(i, 1, len) rev[i] = rev[i >> 1] >> 1 | (i & 1) << (bit - 1);
89
90 F(i, 0, m) {
91 int w = 1LL * a[i] * ifac[i] % P;
92 x[i] = w / p, y[i] = w % p;
93 }
94 F(i, 0, m) {
95 int w = 1LL * b[i] * ifac[i] % P;
96 x2[i] = w / p, y2[i] = w % p;
97 }
98 FFT(x, 0);
99 FFT(y, 0);
100 FFT(x2, 0);
101 FFT(y2, 0);
102 For(i, 0, len) r[i] = x[i] * x2[i];
103 FFT(r, 1);
104 contri(1LL * p * p % P);
105 For(i, 0, len) r[i] = x[i] * y2[i];
106 FFT(r, 1);
107 contri(p);
108 For(i, 0, len) r[i] = y[i] * x2[i];
109 FFT(r, 1);
110 contri(p);
111 For(i, 0, len) r[i] = y[i] * y2[i];
112 FFT(r, 1);
113 contri(1);
114 For(i, 0, len) c[i] = 1LL * c[i] * fac[i] % P;
115
116 int sum = 0;
117 For(i, 0, len) sum = (sum + c[i]) % P;
118 ans = (ans + 1LL * I * sum) % P, printf("%d\n", ans);
119
120 return 0;
121 }
利用递归恒等式展开，可以得到线性做法。
设 $m = n - 2$ 。
设 $\begin{aligned} g_i = \sum_{j = 0} ^ m \binom{i + j}{i} b ^ j \end{aligned}$ 。
贡献为 $\begin{aligned} \sum_{i = 0} ^ m a ^ i g_i \end{aligned}$ 。
递推求 $g$ 。
当 $b = 1$ 时，根据上指标求和恒等式，有 $\begin{aligned} g_i = \sum_{j = 0} ^ m \binom{i + j}{i} = \binom{i + m + 1}{i + 1} \end{aligned}$ 。
当 $b \ne 1$ 时，边界为 $\begin{aligned} g_0 = \sum_{j = 0} ^ m b ^ j \end{aligned}$ ，转移为
$\begin{aligned} g_i & = \sum_{j = 0} ^ m \binom{i + j}{i} b ^ j \\ & = \sum_{j = 0} ^ m \binom{i - 1 + j}{i - 1} b ^ j + \sum_{j = 0} ^ {m - 1} \binom{i + j}{i} b ^ {j + 1} \\ & = g_{i - 1} + b(g_i - \binom{i + m}{i}) b ^ m \\ & = g_{i - 1} + bg_i - b ^ {m + 1} \binom{i + m}{i} \end{aligned}$
解得 $\begin{aligned} g_i = \frac{b ^ {m + 1} \binom{i + m}{i} - g_{i - 1}}{b - 1} \end{aligned}$
实现：
1 #include <bits/stdc++.h>
2 using namespace std;
3 #define F(i, a, b) for (register int i = (a), _b = (b); i <= _b; i++)
4 #define For(i, a, b) for (register int i = (a), _b = (b); i <= _b; i++)
5 inline int rd(void) {
6 int f = 1; char c = getchar(); for (; !isdigit(c); c = getchar()) if (c == ‘-‘) f = -1;
7 int x = 0; for (; isdigit(c); c = getchar()) x = x*10+c-‘0‘; return x*f;
8 }
9
10 const int N = 200005;
11 const int P = 1e6+3;
12
13 int fac[P], inv[P], ifac[P];
14 int n, a[N], b[N], c;
15 int ans;
16 int m, g[N];
17
18 inline int C(int n, int m) {
19 if (m < 0 || n < m) return 0;
20 return 1LL * fac[n] * ifac[m] % P * ifac[n-m] % P;
21 }
22
23 int main(void) {
24 #ifndef ONLINE_JUDGE
25 freopen("4451.in", "r", stdin);
26 #endif
27
28 fac[0] = 1;
29 For(i, 1, P) fac[i] = 1LL * fac[i-1] * i % P;
30 inv[1] = 1;
31 For(i, 2, P) inv[i] = 1LL * (P - P / i) * inv[P % i] % P;
32 ifac[0] = 1;
33 For(i, 1, P) ifac[i] = 1LL * ifac[i-1] * inv[i] % P;
34
35 n = rd(), a[1] = rd(), b[1] = rd(), c = rd();
36 a[0] = b[0] = 1;
37 F(i, 2, n) a[i] = 1LL * a[i-1] * a[1] % P;
38 F(i, 2, n) b[i] = 1LL * b[i-1] * b[1] % P;
39
40 F(i, 1, n) {
41 int w = rd();
42 if (i > 1) ans = (ans + 1LL * w * a[n-1] % P * b[n-i] % P * C(n+n-i-2, n-2)) % P;
43 }
44 F(i, 1, n) {
45 int w = rd();
46 if (i > 1) ans = (ans + 1LL * w * a[n-i] % P * b[n-1] % P * C(n+n-i-2, n-2)) % P;
47 }
48
49 m = n-2;
50 if (b[1] == 1) {
51 F(i, 0, m)
52 g[i] = C(i+m+1, i+1);
53 }
54 else {
55 F(j, 0, m) g[0] = (g[0] + b[j]) % P;
56 F(i, 1, m)
57 g[i] = (1LL * C(i+m, i) * b[m+1] - g[i-1]) % P * inv[b[1] - 1] % P;
58 }
59 int sum = 0;
60 F(i, 0, m) sum = (sum + 1LL * g[i] * a[i]) % P;
61 ans = (ans + 1LL * c * sum) % P, printf("%d\n", ans);
62
63 return 0;
64 }