根号分治 + NTT优化dp转移：ICPC2026网络赛B

image-20260911155646417

考虑一个朴素的dp

定义 cic_i 表示有多少个 aa 小于等于 ii。显然 ii 只有模这 cic_i 个数才会改变。

设 fif_i 表示在整个过程中，有多少种排列使得 ii 存在过。显然对于 i<a1i<a_1，fif_i 就是答案。而 i≥a1i\ge a_1，答案是0。

fi×1ci→j∈{a}fi  mod  jf_i\times \frac 1{c_i}\to_{j\in \{a\}}f_{i\,\bmod\, j} 

1
2
3
4
5
6
f[m] = fac[n]; 
for(v = m / 2; v >= 0; --v) {
    for(i = c[v] + 1; i <= n; ++i) 
        for(u = a[i] + v; u <= m; u += a[i]) 
            Add(f[v], f[u] * inv[c[u]]); 
}

我们发现一个性质。i mod ji\bmod j 只要能成功取模，必然小于 i2\dfrac {i}2。因此这启示我们可以进行分块处理，按 BB 分块后，除了最后一块，其他不可能存在内部之间的转移。

我们考虑现在对于一个块 y∈[l,r]y\in [l,r] 内的转移，我们只考虑从 >r>r 的位置 uu 转移过来。接下来我们考虑模数 x>rx > r。必然满足 u=kx+yu=kx+y，即 x∣u−yx|u-y。令 gi=fici!g_i=\dfrac {f_i}{c_i!}，dbas=#{a[x]∣bas}d_{bas}=\#\{a[x]|bas\}

fy=∑gudu−yf_y=\sum g_{u}d_{u-y}

右边这个东西就是一个卷积的形式，直接NTT即可。这部分的复杂度是 O(Bnlog⁡n)O(Bn\log n)

然后考虑模数 x∈[l,r]x\in [l,r]，在这种情况下直接暴力枚举即可，复杂度 O(B3)O(B^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
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
signed main() {
	#ifdef LOCAL
	  freopen("in.txt", "r", stdin);
	  freopen("out.txt", "w", stdout);
	#endif
//	srand(time(NULL));
//	T = read();
//	while(T--) {
//
//	}
	n = read(); m = read(); int B = 4096; 
	init(max(n, m)); 
	for(i = 1; i <= n; ++i)	a[i] = read(), mp[a[i]] = 1; 
	for(i = j = 0; i <= m; ++i) {
		while(j < n && a[j + 1] <= i) ++j; 
		c[i] = j;  
	} 
	f[m] = fac[n]; g[m] = f[m] * inv[c[m]] % mo; 
	if(a[1] > m) {
		for(i = 0; i <= m; ++i) printf("%lld ", f[i]); 
		return 0; 
	} 
	for(i = 1; i <= n; ++i) if(a[i] <= m) Add(f[m % a[i]], g[m]); 
	r = m / 2; l = (m / 2) / B * B; prex = m; int M = m / 2; 
	while(r >= 0) {
//		debug("[%lld %lld]\n", l, r); 
		for(x = r + 1; x <= prex; ++x) if(mp[x]) {
			for(bas = x; bas <= m; bas += x) d[bas]++; 
		}
		vector<int>G(M - r + 1), D(m - r + 1); 
		for(i = r + 1; i <= M; ++i) G[M - i] = g[i]; 
		for(i = r + 1; i <= m; ++i) D[i - r] = d[i]; 
		auto H = atcoder :: convolution(G, D); 
		for(y = l; y <= r; ++y) if(y + r <= M) Add(f[y], H[M - y - r]); 
		for(y = l; y <= r; ++y) {
			for(i = c[y] + 1; i <= n && a[i] <= r; ++i) {
				for(x = a[i], u = x + y; u <= M; u += x) 
					Add(f[y], g[u]); 
			}
//			if(y > M) debug("%lld %lld\n", y, M); 
			g[y] = f[y] * inv[c[y]] % mo; //G[M - y] = g[y];   
		}
		prex = r; r = l - 1; l -= B; 
	}
	for(v = prex; v >= 0; --v) {
		for(i = c[v] + 1; i <= n && a[i] + v <= prex; ++i) 
			for(u = a[i] + v; u <= prex; u += a[i]) 
				Add(f[v], g[u]); 
		g[v] = f[v] * inv[c[v]] % mo; 
	}
	for(i = 0; i < a[1]; ++i) printf("%lld ", f[i]); 
	for(; i <= m; ++i) printf("0 "); 
	return 0;
}