「BZOJ 4451」「CERC 2015」Frightful Formula

Posted 2020-10-14

　　给 $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 }