二阶差分生成函数优化dp + NTT + cdq优化NTT：CCPC2026网络赛G

Growing Prefix Span - Problem - QOJ.ac

哭了，这道题调了10多个小时，打了无数个补丁，终于调出来了

好，首先一开始肯定是要处理一下 $[L_i,R_i]$，变成合法区间，然后判断一下无解。

然后我们可以列一个很显然的 $O(n^2)$ 的dp。$f(i,j)$ 代表前 $i$ 个数长为 $j$，显然有：

$$f(i,j)=f(i-1,j)\times (j-i+1)+2\sum_{k<j} f_{i-1,k}$$

上面是很自然的，接下来的部分就不太自然了。

我们首先考虑把 $[L,R]$ 的限制抛弃掉，然后我们拿生成函数来表示dp，设：

$$F_i(x)=\sum_{j=0}^{\infty}f(i,j)x^j$$

则：

$$F_i(x)= xF'_{i-1}(x)-(i-1)F_{i-1}(x)+\frac {2x}{1-x}F_{i-1}(x)$$

这里解释一下：

• 我们要构造 $j\cdot f(i,j)$，其中 $j$ 是 $x$ 的指数，所以我们可以用求导把 $j$ 弄出来。然后指数补个1
• $\dfrac{F_i(x)}{1-x}$ 代表求 $F_i(x)$ 的前缀和（可以用泰勒展开或等比数列求和公式来证）。但我们要求的是 $k<j$，而不是 $k\le j$，所以还要再拿个 $x$ 进行平移。

接下来，我们看到式子中有一个 $\dfrac 1{1-x}$，非常丑陋。我们尝试去乘 $(1-x)$ 来使整个式子变干净。

而我们还有一个求导的东西，而求导会降幂。所以我们考虑构造二次项，那样子出来的东西应该是能消掉的。

设：

$$Q_i(x)=\frac{F_i(x)}{(1-x)^2}$$

这里，$F_i(x)$ 就是 $Q_i(x)$ 的二阶差分。如果我们提取项数就有：

$$f_i(j)=(q_i(j)-q_i(j-1))-(q_i(j-1)-q_i(j-2))$$

然后我们代入上面的式子中（$F_i(x)=(1-x)^2Q_i(x)$）：

$$Q_i(x)(1-x)^2=x(-2(1-x)Q_{i}(x)+(1-x)^2Q'_{i-1}(x))-(i-1)(1-x)^2Q_{i-1}(x)+2x(1-x)Q_{i-1}(x)$$

右式中的第一项和第三项消掉，得到：

$$Q_i(x)(1-x)^2=x(1-x)^2Q'_{i-1}(x)-(i-1)(1-x)^2Q_{i-1}(x)$$

同时除去 $(1-x)^2$：

$$Q_i(x)=xQ'_{i-1}(x)-(i-1)Q_{i-1}(x)$$

提取 $[x^j]$：

$$q_{i}(j)=j\cdot q_{i-1}(j)-(i-1)q_{i-1}(j)\\ q_i(j)=(j-i+1)q_{i-1}(j)$$

我们得到了一个非常完美的式子！根据这个式子，我们可以轻易求出 $q_i(j)$，然后求两次差分就可以还原 $f$ 了。

但这有一个问题：如果有 $[L_i,R_i]$ 怎么办？

我们假设我们已经知道前 $i-1$ 的 $q$ 值，我们现在想快速求出 $q_i$。而且前 $i-1$ 的 $q$ 值是满足 $[L,R]$ 的，所以我们只需要考虑当前的 $[L_i,R_i]$。

我们令没有 $[L_i,R_i]$ 时，我们按照上面求出来得到的 $q$ 记为 $q^*_i(x)$。我们试图求出 $q_i(x)$ 和 $q^*_i(x)$ 的关系。

加入限制后，有 \forall j R_i\Rightarrow f_i(j)=0，而剩余的 $f_i(j)$ 不变，我们看回这个式子，思考会对 $q_i(j)$ 有什么影响：

$$f_i(j)=(q_i(j)-q_i(j-1))-(q_i(j-1)-q_i(j-2))$$

显然，根据定义，

$$q_i(j<L_i)=0$$

对于 $L_i\le j\le R_i$，我们可以手动模拟一下，对 $\Delta f$ 求两次积分（前缀和）：

$q_i(j)$ 相比 $q^*_i(j)$ 只差了一个固定的一次函数，也就是：

$$\begin{align*} q_i(L_i\le j\le R_i)&=q^*_i(j)+A_i+jB_i\\[6px] &=(j-i+1)q_{i-1}(j)+A_i+jB_i \end{align*}$$

这个关系式在 $L_i-1$ 和 $L_i-2$ 时依然成立。因此我们可以直接用 $q^*_i(L_i-1)+A_i+(L_i-1)B_i=0$ 和 $q^*_i(L_i-2)+A_i+(L_i-2)B_i=0$ 即可解出 $A_i$ 和 $B_i$。

考虑 $q_i(j>R_i)$ 是什么情况。这个时候其二阶差分为0，一个多项式的二阶导为0那么它只能是一次函数。因此 $q_i(j>R_i)$ 就是一个一次函数，我们可以直接表示为：

$$q_i(j>R_i)=C_i+jD_i$$

由于在 $R_i$ 和 $R_i-1$ 时也满足此式子，而 $A_i,B_i$ 已知，我们就可以求出 $q_i(R_i)$ 和 $q_i(R_i-1)$，然后通过解方程就可以求出 $C_i,D_i$ 了。

接下里我们回顾一下，哪些 $q_i(j)$ 是我们真正关心的。

不难发现，在上述过程中，对于每个 $i$，其实我们只需要求4个 $j$，分别是：

$$q_{i-1}(L_i-1),q_{i-1}(L_i-2),q_{i}(R_i-1),q_i(R_i)$$

因此我们不需要每个 $q_i(j)$ 都求出来，我们只需要求出我们真正关心的 $q_i(j)$ 即可，因此我们可以考虑求 $q_i(j)$ 的通项。

显然 jR_i 是简单的，我们考虑 $L_i\le j\le R_i$ 的情况。

接下来我们看回这个式子：

$$q_i(L_i\le j\le R_i)=(j-i+1)q_{i-1}(j)+A_i+jB_i$$

我们设法直接求其通项。

不妨令 $e_j$ 表示第一个满足 $j\le R_i$ 的 $i$，其中 $e_j$ 是可以通过双指针非常轻松预处理的。直接展开递归式，我们有：

$$q_i(j)=\sum_{k=e_j}^i(A_k+jB_k)\prod_{t=k+1}^i(j-t+1)+q_{e_j-1}(j)\prod_{t=e_j}^i(j-t+1)$$

令 $H_j=q_{e_j-1}(j)$，则 $H_j=C_{e_j-1}+jD_{e_j-1}$，这个我们在 $i=e_j-1$ 的时候就可以直接求出来的。然后代换进去，同时用阶乘处理连乘：

$$q_i(j)=\sum_{k=e_j}^i(A_k+jB_k)\cdot \frac{(j-k)!}{(j-i)!}+H(j)\frac{(j-e_j+1)!}{(j-i)!}\\[6px] q_i(j)=\frac1{(j-i)!}\left(\sum_{k=e_j}^i(A_k+jB_k)\cdot (j-k)!+H(j)(j-e_j+1)!\right)$$

也就是说，只要我们知道了 $A,B,C,D,H$，我们就可以求出任意的 $q_i(j)$。根据上面的分析，这5个变量都是可以求的。我们可以轻松的写出一份达到 $O(n^3)$ 的暴力：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26

for(i = j = 1; i <= n; ++i) {
    while(j <= n && R[j] < i) ++j; 
    e[i] = j; 
}
for(j = 1; j <= n; ++j) q[1][j] = j; 
for(j = R[1] + 1; j <= R[2]; ++j) H[j] = q[1][j]; 
for(i = 2; i <= n; ++i) {
    u = ((L[i] - 1) - i + 1) * q[i - 1][L[i] - 1]; Mod(u); 
    v = ((L[i] - 2) - i + 1) * q[i - 1][L[i] - 2]; Mod(v); 
    B[i] = u - v; A[i] = u - B[i] * (L[i] - 1); 
    A[i] = -A[i]; B[i] = -B[i]; 
    Mod(A[i]); Mod(B[i]); 
    for(j = L[i]; j <= R[i]; ++j) {
        q[i][j] = H[j] * fac[j - e[j] + 1]; 
        for(k = e[j]; k <= i; ++k) 
            Add(q[i][j], (A[k] + B[k] * j) % mo * fac[j - k]); 
        Mul(q[i][j], ifac[j - i]); 
    }
    u = q[i][R[i]]; v = q[i][R[i] - 1]; 
    D[i] = u - v; C[i] = u - D[i] * R[i]; 
    for(j = R[i] + 1; j <= n; ++j)
        q[i][j] = C[i] + D[i] * j, Mod(q[i][j]); 
    for(j = R[i] + 1; j <= R[i + 1]; ++j) H[j] = q[i][j]; 
}
ans = C[n] + D[n] * n; Mod(ans); 
printf("%lld", ans); 

所以剩下的问题是：在已知 $A,B,C,D,H$ 的情况下，如何加速求 $q_i(j)$ 的过程？

$(j-i)!$ 直接扔掉。$H(j)(j-e_j+1)!$ 可以O(1)计算，也扔掉，然后我们剩余的本质上在求这个：

$$\sum_{k=e_j}^i(A_k+jB_k)\cdot (j-k)!\\ =\sum_{k=e_j}^iA_k\cdot (j-k)!+j\cdot \sum_{k=e_j}^i B_k\cdot (j-k)!$$

此时 $A_k$ 和 $B_k$ 的求法等价了。现在考虑具体怎么求其中任意一个。

这个东西我们可以写成卷积的形式，令 $F_j(x)=\sum_{k=e_j}^{i}A_kx^k,G_j(x)=\sum_{k=e_j}^i(j-k)!x^{j-k}$，然后我们就相当于求 $[x^j]F_j\times G_j$

利用卷积，我们希望做到的是，能一口气求出很多个 $j$ 的答案。但遗憾的是，上面这个卷积有 $[e_j,i]$ 的限制。因此我们可以考虑采用分段卷积。

考虑cdq分治。首先因为有一个分治的存在，所以 $i$ 的上界就可以直接扔了。等我们求到 $i$ 的时候直接读取当前的 $[x^j]F_j\times G_j$ 就是答案。（这里有一个细节是求 $q_{i-1}(L_i-1),q_{i-1}(L_i-2)$ 要在 $i-1$ 的时候提前求好，不能等到 $i$ 的时候再求）

但这个下界 $e_j$ 有点恶心，接下来我介绍我的方法来保证所有贡献不重不漏：

1. 不重

   考虑我们现在已经完成了 $[l,r]$ 整个区间的cdq了。在这个时候，我对 $[mid+1,r]$ 内的所有数进行卷积贡献。而我只统计满足 $e_j\in[l,mid]$ 的 $j$，也就是只和这些 $j$ 进行上面的卷积。这样子显然是每一层右区间只会和当前层左区间的合法 $j$ 转移，因此是不重的。

   同时，由于 $e_j$ 是单增的，首先合法的 $j$ 必然是连续一段区间，因此这个区间可以直接拿二分求出。更进一步，那样子合法 $j$ 的长度均摊下来复杂度就是对的了。$[mid+1,r]$ 均摊长度是 $\log n$，$j$ 也是 $\log n$。因此卷起来就是 $\log^2n$。（分治每一层所有NTT的长度加起来严格为 $O(n)$）

2. 不漏

   但我们上面这个做法有一个问题。我们现在考虑我们处理一个区间 $[l,r]$，左区间 $[l,mid]$ 对满足 $e_j<l$ 的 $j$ 的贡献没有给右区间。而这些 $j$ 又太多了，我们显然不能全部弄，那怎么办呢？

   我们发现，因为 $L_i,R_i$ 是单增的（预处理后）。因此我们在 $[mid+1,r]$ 时我们所询问的 $j$ 必然在 $j\in[L(mid+2)-2,L(r+1)-1]\cup[R(mid+1)-1,R(r)]$ 中，因此我们真正需要的就是这两个区间里的 $j$ （具体还要处理一下交并的细节）。所以我们卷积的 $j$ 通过右区间的需求来限制住就好。当然，也必须满足 $e_j<l$。

   此时我们是拿右边的区间来限制，所以均摊下来每层NTT的长度也是严格 $O(n)$ 的。

   需要注意的是，这一部分的贡献在处理完右区间之后是要重新减掉的。

由于总共 $\log n$ 层，每层NTT的长度严格 $O(n)$，因此总复杂度为 $O(n\log^2n)$.

题解讲完q的定义后截断、cdq和NTT的部分就没看懂了，然后自己手推了一大轮，慢慢从暴力往上一步步改，边对拍边打补丁，力竭了…

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162

using namespace std;bool Mbe;
using namespace atcoder; 
#ifdef LOCAL
 #define debug(...) fprintf(stdout, ##__VA_ARGS__)
#else
 #define debug(...) void(0)
#endif
#define int long long
inline int read(){int x=0,f=1;char ch=getchar();
while(ch<'0'||ch>'9'){if(ch=='-')f=-1;
ch=getchar();}while(ch>='0'&&ch<='9'){x=(x<<1)+
(x<<3)+(ch^48);ch=getchar();}return x*f;}
#define Z(x) (x)*(x)
#define pb push_back
#define fi first
#define se second
#define mo 998244353
//#define M
#define N 200010
int pw(int a, int b) {
	int ans = 1; 
	while(b) {
		if(b & 1) ans *= a; 
		a *= a; b >>= 1; 
		ans %= mo; a %= mo; 
	}
	return ans; 
}
int pw(int a) { return pw(a, mo - 2); }
int fac[N], inv[N], ifac[N]; 
void init(int n) {
	int i; 
	for(i = fac[0] = 1; i <= n; ++i) fac[i] = fac[i-1] * i % mo; 
	ifac[n] = pw(fac[n], mo-2); 
	for(i = n - 1; i >= 0; --i) ifac[i] = ifac[i+1] * (i+1) % mo; 
    for(i = 1; i <= n; ++i) inv[i] = ifac[i] * fac[i-1] % mo; 
}
inline void Mod(int &a) { if(a >= mo || a <= -mo) a %= mo; if(a < 0) a += mo; }
inline void Add(int &a, int b) { a += b; Mod(a); }
inline void Mul(int &a, int b) { Mod(b); a *= b; Mod(a); } 
const int iv2 = pw(2, mo - 2); 
int n, m, i, j, k, T;
int L[N], R[N]; 
int A[N], B[N], C[N], D[N], e[N], H[N], u, v, ans; 
int fa[N], fb[N], FL1[N], FL2[N]; 
//vector<int>V[N];  

int getQ(int i, int j) {
	if(i == 1) return j; 
	if(j < L[i]) return 0; 
	if(j > R[i]) return (C[i] + D[i] * j) % mo; 
	int Q = H[j] * fac[j - e[j] + 1] % mo, idx; 
	Add(Q, fa[j] + fb[j] * j); 
	Mul(Q, ifac[j - i]); 
	Add(Q, A[i] + B[i] * j); 
	return Q; 
}

void work(int l, int mid, int r) {
	vector<int>a, b, f; 
	int j1 = lower_bound(e + 1, e + n + 1, l + 1) - e; 
	int j2 = upper_bound(e + 1, e + n + 1, mid + 1) - e - 1; 
	if(j1 > j2 || j1 > n || e[j1] - 1 > mid) return ; 
	for(k = mid + 1; k <= r; ++k) a.pb(A[k]), b.pb(B[k]); 
	for(k = j1 - r; k <= j2 - (mid + 1); ++k) 
		if(k >= 0) f.pb(fac[k]); 
		else f.pb(0); 
	a = convolution(a, f); 
	b = convolution(b, f); 
	for(j = j1, i = r - (mid + 1); j <= j2; ++j, ++i) {
		Add(fa[j], a[i]), 
		Add(fb[j], b[i]); 
	}
		
}

void contribute(int l, int mid, int L, int R, int o) {
	vector<int>a, b, f; 
	if(L == 0 || R == 0 || L > R) return ; 
	int j1 = lower_bound(e + 1, e + n + 1, 1 + 1) - e; 
	int j2 = upper_bound(e + 1, e + n + 1, l - 1 + 1) - e - 1; 
	j1 = max(j1, L); j2 = min(j2, R); 
	if(j1 > j2 || j1 > n) return ; 
	for(k = l; k <= mid; ++k) a.pb(A[k]), b.pb(B[k]); 
	for(k = j1 - mid; k <= j2 - l; ++k) 
		if(k >= 0) f.pb(fac[k]); 
		else f.pb(0); 
	a = convolution(a, f); 
	b = convolution(b, f); 
	for(j = j1, i = mid - l; j <= j2; ++j, ++i) {
		Add(fa[j], a[i] * o), 
		Add(fb[j], b[i] * o); 
	}
}


pair<pair<int, int>, pair<int, int> > Devide(int l1, int r1, int l2, int r2) {
	if(l2 <= 0 || r2 <= 0) return {{l1, r1}, {0, 0}}; 
	if(l2 <= r1 && l1 <= r2) return {{min(l1, l2), max(r1, r2)}, {0, 0}}; 
	return {{l1, r1}, {l2, r2}}; 
}

void solve(int l, int r) {
	if(l == r) {
		if(l == 1) {
			FL1[1] = getQ(1, L[2] - 1); FL2[1] = getQ(1, L[2] - 2); 
			return ; 
		}
		int i = l; 
		u = ((L[i] - 1) - i + 1) * FL1[i - 1]; Mod(u); 
		v = ((L[i] - 2) - i + 1) * FL2[i - 1]; Mod(v); 
		B[i] = u - v; A[i] = u - B[i] * (L[i] - 1); 
		A[i] = -A[i]; B[i] = -B[i]; 
		Mod(A[i]); Mod(B[i]); 
		u = getQ(i, R[i]); v = getQ(i, R[i] - 1); 
		D[i] = u - v; C[i] = u - D[i] * R[i]; 
		for(j = R[i] + 1; j <= R[i + 1]; ++j) H[j] = getQ(i, j); 
		if(i < n) FL1[i] = getQ(i, L[i + 1] - 1), FL2[i] = getQ(i, L[i + 1] - 2); 
		return ; 
	}
	int mid = (l + r) >> 1; 
	solve(l, mid); 
	auto [t1, t2] = Devide(R[mid + 1] - 1, R[r], L[mid + 2] - 2, L[r + 1] - 1); 
	contribute(l, mid, t1.first, t1.second, 1); 
	contribute(l, mid, t2.first, t2.second, 1); 
	solve(mid + 1, r); 
	contribute(l, mid, t1.first, t1.second, -1); 
	contribute(l, mid, t2.first, t2.second, -1); 
	work(l, mid, r); 
}

bool Med; signed main() {
	#ifdef LOCAL
	  freopen("in.txt", "r", stdin);
	  freopen("out.txt", "w", stdout);
	  atexit([]{fprintf(stderr,"%.lfms\n%lfMB\n",clock()*1000./CLOCKS_PER_SEC,(&Mbe-&Med)/1024./1024);});
	#endif
//	srand(time(NULL));
//	T = read();
//	while(T--) {
//
//	}
	n = read(); init(n); 
	for(i = 1; i <= n; ++i) L[i] = read(), R[i] = read(); 
	for(i = 1; i <= n; ++i)
		L[i] = max({L[i], L[i - 1], i});  
	for(i = n - 1; i >= 1; --i) R[i] = min(R[i], R[i + 1]); 
	R[1] = 1; 
	for(i = 1; i <= n; ++i)
		if(L[i] > R[i]) { printf("0"); return 0; }
	if(n == 1) { printf("1"); return 0; }
	for(i = j = 1; i <= n; ++i) {
		while(j <= n && R[j] < i) ++j; 
		e[i] = j; 
	}
	L[n + 1] = R[n + 1] = n + 1; 
	for(j = R[1] + 1; j <= R[2]; ++j) H[j] = getQ(1, j); 
	solve(1, n); 
	ans = C[n] + D[n] * n; Mod(ans); 
	printf("%lld", ans); 
	return 0;
}