[EZOI][1211NOI模拟赛]LCM(根号分治+容斥)

§ 1 题意

给定一个大小为 $n$ 的数集 $S$,求 $\sum\limits_{S’ \subseteq S,S’ \neq \emptyset} {\rm LCM}(S’)$,即 $S$ 所有非空子集的 ${\rm LCM}$ 和。

$n \leq 2000,a_i \leq 200$。


§ 2 分析

首先发现 $a_i \leq 200$ 很小,考虑是否可以分解质因数统计。但 $200$ 以内的质数有 $46$ 个,组成的 ${\rm LCM}$ 值域过大。

考虑根号分治,观察发现 $\gt \sqrt{a_i} \approx 14$ 的质数的幂次最高为 $1$,而 $\leq 14$ 的质数只有 $6$ 个,分别为 $2,3,5,7,11,13$,其最高幂次分别为 $7,4,3,2,2,2$。

所以若只考虑 $2,3,5,7,11,13$,则 ${\rm LCM}$ 的不同值至多只有 $8\times5\times4\times3\times3\times3=4320$ 种。

我们将 $a_i$ 按除完 $2,3,5,7,11,13$ 后剩下的值 $c_i$ 分类,可以证明 $c_i=1 $ 或 $\gt 13$ 的质数。

令 $f_{c,p_1,p_2,p_3,p_4,p_5,p_6}$ 表示除完后剩下的值为 $c$ 时,是 $2^{p_1}\times3^{p_2}\times5^{p_3}\times7^{p_4}\times11^{p_5}\times13^{p_6}$ 因数的数的个数:

1
2
3
4
5
6
7
8
for(register int i = 1; i <= n; i++)
for(register int p1 = a[i][1]; p1 <= mx[1]; p1++)
for(register int p2 = a[i][2]; p2 <= mx[2]; p2++)
for(register int p3 = a[i][3]; p3 <= mx[3]; p3++)
for(register int p4 = a[i][4]; p4 <= mx[4]; p4++)
for(register int p5 = a[i][5]; p5 <= mx[5]; p5++)
for(register int p6 = a[i][6]; p6 <= mx[6]; p6++)
f[c[i]][p1][p2][p3][p4][p5][p6]++;

令 $g_{p1,p2,p3,p4,p5,p6}$ 表示不考虑 $\gt 13$ 的因子时 ${\rm LCM}$ 为 $2^{p_1}\times3^{p_2}\times5^{p_3}\times7^{p_4}\times11^{p_5}\times13^{p_6}$ 的因数的子集数。

首先考虑除完后剩下的值 $c=1$ 时,$F=f_{1,p_1,p_2,p_3,p_4,p_5,p_6}$ 个数中的每个要么选要么不选,共有 $2^F$ 种可能:

1
2
3
4
5
6
7
for(register int p1 = 0; p1 <= mx[1]; p1++)
for(register int p2 = 0; p2 <= mx[2]; p2++)
for(register int p3 = 0; p3 <= mx[3]; p3++)
for(register int p4 = 0; p4 <= mx[4]; p4++)
for(register int p5 = 0; p5 <= mx[5]; p5++)
for(register int p6 = 0; p6 <= mx[6]; p6++)
g[p1][p2][p3][p4][p5][p6] = bin[f[1][p1][p2][p3][p4][p5][p6]];

然后考虑除完后剩下的值 $c\gt 13$ 时,$F=f_{c,p_1,p_2,p_3,p_4,p_5,p_6}$ 个数中,都不选有 $1$ 种可能,选至少一个都会导致 ${\rm LCM}$ 扩大 $c$ 倍,由于只需统计 ${\rm LCM}$ 的和,可将一种算作 $c$ 种,等价于共有 $1+c\times(2^F-1)$ 种可能:

1
2
3
4
5
6
7
8
9
10
for(register int i = 1; i <= 40; i++){
const int &r = rem[i];
for(register int p1 = 0; p1 <= mx[1]; p1++)
for(register int p2 = 0; p2 <= mx[2]; p2++)
for(register int p3 = 0; p3 <= mx[3]; p3++)
for(register int p4 = 0; p4 <= mx[4]; p4++)
for(register int p5 = 0; p5 <= mx[5]; p5++)
for(register int p6 = 0; p6 <= mx[6]; p6++)
(g[p1][p2][p3][p4][p5][p6] *= (r * (bin[f[r][p1][p2][p3][p4][p5][p6]] - 1LL + MOD) + 1LL) % MOD) %= MOD;
}

接下来通过容斥原理即可将 $g_{p1,p2,p3,p4,p5,p6}$ 转化为不考虑 $\gt 13$ 的因子时 ${\rm LCM}$ 严格等于 $2^{p_1}\times3^{p_2}\times5^{p_3}\times7^{p_4}\times11^{p_5}\times13^{p_6}$ 的子集数,注意容斥要从大到小进行:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
for(register int p1 = mx[1]; ~p1; p1--)
for(register int p2 = mx[2]; ~p2; p2--)
for(register int p3 = mx[3]; ~p3; p3--)
for(register int p4 = mx[4]; ~p4; p4--)
for(register int p5 = mx[5]; ~p5; p5--)
for(register int p6 = mx[6]; ~p6; p6--)
for(register int d1 = 0; d1 <= 1 && d1 <= p1; d1++)
for(register int d2 = 0; d2 <= 1 && d2 <= p2; d2++)
for(register int d3 = 0; d3 <= 1 && d3 <= p3; d3++)
for(register int d4 = 0; d4 <= 1 && d4 <= p4; d4++)
for(register int d5 = 0; d5 <= 1 && d5 <= p5; d5++)
for(register int d6 = 0; d6 <= 1 && d6 <= p6; d6++) if(d1 | d2 | d3 | d4 | d5 | d6){
if(d1 ^ d2 ^ d3 ^ d4 ^ d5 ^ d6)
(g[p1][p2][p3][p4][p5][p6] += MOD - g[p1 - d1][p2 - d2][p3 - d3][p4 - d4][p5 - d5][p6 - d6]) %= MOD;
else (g[p1][p2][p3][p4][p5][p6] += g[p1 - d1][p2 - d2][p3 - d3][p4 - d4][p5 - d5][p6 - d6]) %= MOD;
}

最后答案即为 $\sum\limits_{p1,p2,p3,p4,p5,p6}2^{p_1}\times3^{p_2}\times5^{p_3}\times7^{p_4}\times11^{p_5}\times13^{p_6}\times g_{p1,p2,p3,p4,p5,p6}$ 减去 $1$,即除去空集。

总时间复杂度为 $O(8\times5\times4\times3\times3\times3\times(n+40+2^6))​$。


§ 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
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
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<cmath>
#include<iostream>
#include<algorithm>
#define MOD 1000000007
using namespace std;
typedef long long ll;
const int prm[7] = {1, 2, 3, 5, 7, 11, 13};
const int rem[41] = {1, 17, 19, 23, 29, 31, 37, 41, 43, 47, 53,
59, 61, 67, 71, 73, 79, 83, 89, 97, 101,
103, 107, 109, 113, 127, 131, 137, 139, 149, 151,
157, 163, 167, 173, 179, 181, 191, 193, 197, 199};

inline void getint(int &num){
register int ch, neg = 0;
while(!isdigit(ch = getchar())) if(ch == '-') neg = 1;
num = ch & 15;
while(isdigit(ch = getchar())) num = num * 10 + (ch & 15);
if(neg) num = -num;
}

int n, mx[7], a[2005][7], c[2005], v[8][5][4][3][3][3], f[205][8][5][4][3][3][3], bin[2005], ans;
ll g[8][5][4][3][3][3];

int main(){
getint(n), bin[0] = 1;
for(register int i = 1; i <= n; i++) bin[i] = (bin[i - 1] << 1) % MOD;
for(register int i = 1; i <= n; i++){
int v; getint(v);
for(register int j = 1; j <= 6; j++)
while(v % prm[j] == 0){
if(++a[i][j] > mx[j]) mx[j] = a[i][j];
v /= prm[j];
}
c[i] = v;
}
for(register int p1 = 0; p1 <= mx[1]; p1++)
for(register int p2 = 0; p2 <= mx[2]; p2++)
for(register int p3 = 0; p3 <= mx[3]; p3++)
for(register int p4 = 0; p4 <= mx[4]; p4++)
for(register int p5 = 0; p5 <= mx[5]; p5++)
for(register int p6 = 0; p6 <= mx[6]; p6++){
if(p1) v[p1][p2][p3][p4][p5][p6] = v[p1 - 1][p2][p3][p4][p5][p6] * 2 % MOD;
else if(p2) v[p1][p2][p3][p4][p5][p6] = v[p1][p2 - 1][p3][p4][p5][p6] * 3 % MOD;
else if(p3) v[p1][p2][p3][p4][p5][p6] = v[p1][p2][p3 - 1][p4][p5][p6] * 5 % MOD;
else if(p4) v[p1][p2][p3][p4][p5][p6] = v[p1][p2][p3][p4 - 1][p5][p6] * 7 % MOD;
else if(p5) v[p1][p2][p3][p4][p5][p6] = v[p1][p2][p3][p4][p5 - 1][p6] * 11 % MOD;
else if(p6) v[p1][p2][p3][p4][p5][p6] = v[p1][p2][p3][p4][p5][p6 - 1] * 13 % MOD;
else v[p1][p2][p3][p4][p5][p6] = 1;
} // 预处理 2^p1 * 3^p2 * 5^p3 * 7^p4 * 11^p5 * 13^p6 的值
for(register int i = 1; i <= n; i++)
for(register int p1 = a[i][1]; p1 <= mx[1]; p1++)
for(register int p2 = a[i][2]; p2 <= mx[2]; p2++)
for(register int p3 = a[i][3]; p3 <= mx[3]; p3++)
for(register int p4 = a[i][4]; p4 <= mx[4]; p4++)
for(register int p5 = a[i][5]; p5 <= mx[5]; p5++)
for(register int p6 = a[i][6]; p6 <= mx[6]; p6++)
f[c[i]][p1][p2][p3][p4][p5][p6]++;
// 按 c[i] 分组统计可能对 LCM = 2^p1 * 3^p2 * 5^p3 * 7^p4 * 11^p5 * 13^p6 有贡献的数的个数
for(register int p1 = 0; p1 <= mx[1]; p1++)
for(register int p2 = 0; p2 <= mx[2]; p2++)
for(register int p3 = 0; p3 <= mx[3]; p3++)
for(register int p4 = 0; p4 <= mx[4]; p4++)
for(register int p5 = 0; p5 <= mx[5]; p5++)
for(register int p6 = 0; p6 <= mx[6]; p6++)
g[p1][p2][p3][p4][p5][p6] = bin[f[1][p1][p2][p3][p4][p5][p6]];
// c[i] = 1 的对答案贡献倍数为 2^f
for(register int i = 1; i <= 40; i++){
const int &r = rem[i];
for(register int p1 = 0; p1 <= mx[1]; p1++)
for(register int p2 = 0; p2 <= mx[2]; p2++)
for(register int p3 = 0; p3 <= mx[3]; p3++)
for(register int p4 = 0; p4 <= mx[4]; p4++)
for(register int p5 = 0; p5 <= mx[5]; p5++)
for(register int p6 = 0; p6 <= mx[6]; p6++)
(g[p1][p2][p3][p4][p5][p6] *= (r * (bin[f[r][p1][p2][p3][p4][p5][p6]] - 1LL + MOD) + 1LL) % MOD) %= MOD;
} // c[i] > 13 的对答案贡献倍数为 1 + (c[i] * 2^f - 1)
for(register int p1 = mx[1]; ~p1; p1--)
for(register int p2 = mx[2]; ~p2; p2--)
for(register int p3 = mx[3]; ~p3; p3--)
for(register int p4 = mx[4]; ~p4; p4--)
for(register int p5 = mx[5]; ~p5; p5--)
for(register int p6 = mx[6]; ~p6; p6--)
for(register int d1 = 0; d1 <= 1 && d1 <= p1; d1++)
for(register int d2 = 0; d2 <= 1 && d2 <= p2; d2++)
for(register int d3 = 0; d3 <= 1 && d3 <= p3; d3++)
for(register int d4 = 0; d4 <= 1 && d4 <= p4; d4++)
for(register int d5 = 0; d5 <= 1 && d5 <= p5; d5++)
for(register int d6 = 0; d6 <= 1 && d6 <= p6; d6++) if(d1 | d2 | d3 | d4 | d5 | d6){
if(d1 ^ d2 ^ d3 ^ d4 ^ d5 ^ d6)
(g[p1][p2][p3][p4][p5][p6] += MOD - g[p1 - d1][p2 - d2][p3 - d3][p4 - d4][p5 - d5][p6 - d6]) %= MOD;
else (g[p1][p2][p3][p4][p5][p6] += g[p1 - d1][p2 - d2][p3 - d3][p4 - d4][p5 - d5][p6 - d6]) %= MOD;
} // 通过容斥将"小于等于"转化为严格等于
for(register int p1 = 0; p1 <= mx[1]; p1++)
for(register int p2 = 0; p2 <= mx[2]; p2++)
for(register int p3 = 0; p3 <= mx[3]; p3++)
for(register int p4 = 0; p4 <= mx[4]; p4++)
for(register int p5 = 0; p5 <= mx[5]; p5++)
for(register int p6 = 0; p6 <= mx[6]; p6++)
ans = (ans + g[p1][p2][p3][p4][p5][p6] * v[p1][p2][p3][p4][p5][p6]) % MOD;
printf("%d\n", (ans + MOD - 1) % MOD);
return 0;
}