题目大意

给出多重数集 S=x1,x2,⋯ ,xnS={x_1,x_2,\cdots,x_n}，求其所有大小为 kk 的子集的 gcd⁡\gcd 之积。共 TT 组数据，要求答案  mod P\bmod P 输出。

不保证 SS 中所有数各不相同，不保证 PP 为质数。T≤60,106≤P≤1014,1≤xi≤8×104,1≤∣S∣≤40000,1≤k≤min⁡(∣S∣,30)T≤60,10^6≤P≤10^{14},1≤x_i≤8\times 10^4,1≤|S|≤40000,1≤k≤\min(|S|,30)。

思路

枚举所有可能出现的 gcd⁡\gcd，设 gcd⁡=g\gcd=g 的长度为 kk 的子集个数为 cgc_g，那么最终答案 $ans=\prod\limits_{g=1}^{\max\limits_{1\le i\le n}\{x_i\}}g^{c_g}$。我们的任务是求出每一个 cgc_g。

枚举 gg 时，统计 SS 中 gg 的倍数出现的个数 numg=∑i=1n[g∣xi]num_g=\sum\limits_{i=1}^n[g\mid x_i]，从这 numgnum_g 个中数字中选 kk 个组成备选集合 S′S'，共 CnumgkC_{num_g}^k 中方案，S′S' 中所有数的 gcd⁡\gcd 必然是 gg 的倍数，从 CnumgkC_{num_g}^{k} 个方案中去掉 gcd⁡=2g,3g⋯\gcd=2g,3g\cdots 的方案，剩下的方案数即为 cgc_g。只要倒叙枚举 gg，cg=Cnumgk−c2g−c3g−⋯c_g=C_{num_g}^k-c_{2g}-c_{3g}-\cdots。

题目要求答案对 PP 取模，考虑到 cgc_g 很大，我们要进行欧拉降幂。需要注意的是，当 gcd⁡(a,m)=1\gcd(a,m)=1，也可直接用下面式子的第三条。

$$ 扩展欧拉定理\quad a^b \equiv \begin{cases} a^{b\ \bmod\ \varphi(m)}, &\gcd(a,m) = 1, \\ a^b, &\gcd(a,m)\ne 1, b < \varphi(m),\\ a^{b\ \bmod\ \varphi(m) +\varphi(m)}, &\gcd(a,m)\ne 1, b \ge \varphi(m). \end{cases} \pmod m $$
$$ 欧拉函数\quad\varphi(m)=m\times \prod\limits_{p_i\in \mathbb{prime},p_i\mid m}\left(1-\frac{1}{p_i}\right) $$

考虑到 PP 很大，我们用 Pollard-Rho 算法质因数分解得到 PP 的所有质因子，再计算 φ(P)\varphi(P)。同时由于最后统计答案多次用到组合数，且组合数是作为指数出现，因此可以预处理组合数；因为 φ(P)φ (P) 不一定为质数，同时 kk 不大，所以我们可以用杨辉三角递推组合数。模乘法需要用龟速乘，但是有__int128。

最终要枚举 g∈[1,max⁡{xi}]g\in[1,\max\{x_i\}] 计算答案，对当前于 gg，要枚举 c2g,c3g,⋯c_{2g},c_{3g},\cdots 计算 cgc_g，总的时间复杂度是 $O\left(\max\limits_{1\le i\le n}\{x_i\}\ln\left(\max\limits_{1\le i\le n}\{x_i\}\right)\right)$。

代码

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
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
#include<iostream>
#include<ctime>
#include<algorithm>
#include<cstring>
#include<unordered_set>
using namespace std;
typedef long long ll;
typedef __int128 bll;
int _;
unordered_set<ll>primeFactor;
int n, k;
ll P, phiP;
ll combine[80004][32];
int cnt[80004];
ll c[80004];
// 快速幂
ll QuickPow(ll a, ll b, ll MOD)
{
	a %= MOD;
	ll res = 1 % MOD;
	while (b)
	{
		if (b & 1)
		{
			res = (bll)res * (bll)a % MOD;
		}
		b >>= 1;
		a = (bll)a * (bll)a % MOD;
	}
	return res;
}
// 素数测试
bool MillerRabin(ll n, int repeat)
{
	if (n == 2 || n == 3)
	{
		return true;
	}
	else if (n % 2 == 0 || n == 1)
	{
		return false;
	}
	// 将n-1分解成2^s*d
	ll d = n - 1;
	int s = 0;
	while (!(d & 1))
	{
		++s;
		d >>= 1;
	}
	while (repeat--)
	{
		// 取随机数[2,n-2]
		ll a = rand() % (n - 3) + 2;
		ll x = QuickPow(a, d, n);
		ll y = 0;
		for (int i = 1;i <= s;i++)
		{
			y = (bll)x * (bll)x % n;
			if (y == 1 && x != 1 && x != n - 1)
			{
				return false;
			}
			x = y;
		}
		if (y != 1)
		{
			return false;
		}
	}
	return true;
}
// 返回最大质因子
ll PollardRho(ll n, ll c)
{
	ll x = rand() % (n - 2) + 1;
	ll y = x, i = 1, k = 2;
	while (1)
	{
		++i;
		x = (bll)x * (bll)x % n + c + n;
		ll d = abs(__gcd(y - x, n));
		if (1 < d && d < n)
		{
			return d;
		}
		else if (y == x)
		{
			return n;
		}
		else if (i == k)
		{
			y = x;
			k <<= 1;
		}
	}
}
// 分解质因数
void FindPrimeFactor(ll n, ll c)
{
	if (n == 1)
	{
		return;
	}
	if (MillerRabin(n, 50))
	{
		primeFactor.emplace(n);
		return;
	}
	ll p = n;
	while (p >= n)
	{
		p = PollardRho(p, c--);
	}
	FindPrimeFactor(p, c);
	FindPrimeFactor(n / p, c);
}
void init()
{
	// 求P的欧拉函数值
	FindPrimeFactor(P, 103);
	phiP = P;
	for (const ll& p : primeFactor)
	{
		phiP = phiP / p * (p - 1);
	}

	// 预处理组合数
	for (int i = 1;i <= n;i++)
	{
		combine[i][0] = combine[i][i] = 1;
		for (int j = 1;j < i && j <= k;j++)
		{
			// 杨辉三角递推
			combine[i][j] = combine[i - 1][j] + combine[i - 1][j - 1];
			// 欧拉降幂
			if (combine[i][j] >= phiP)
			{
				combine[i][j] = combine[i][j] % phiP + phiP;
			}
		}
	}
}
int main()
{
	srand(time(0));
	for (cin >> _;_;_--)
	{
		scanf("%d %d %lld", &n, &k, &P);
		init();
		int maxX = 0;
		for (int i = 1;i <= n;i++)
		{
			int x;
			scanf("%d", &x);
			maxX = max(maxX, x);
			// 预处理cnt数组 记录x出现的次数
			++cnt[x];
		}
		ll ans = 1;
		for (int g = maxX;g >= 1;g--)
		{
			int num = 0;
			// 得到g的倍数的个数
			for (int x = g;x <= maxX;x += g)
			{
				num += cnt[x];
			}
			c[g] = combine[num][k];
			// 排除公约数为g的倍数的集合
			for (int x = 2 * g;x <= maxX;x += g)
			{
				c[g] -= c[x];
			}
			if (c[g] >= phiP)
			{
				c[g] = c[g] % phiP + phiP;
			}
			else if (c[g] < 0)
			{
				c[g] = (c[g] % phiP + phiP) % phiP;
			}
			ans = (bll)ans * (bll)QuickPow(g, c[g], P) % P;
		}
		printf("%lld\n", ans);
		memset(cnt, 0, sizeof(int) * (maxX + 1));
		memset(c, 0, sizeof(ll) * (maxX + 1));
		primeFactor.clear();
	}
	return 0;
}