2021牛客暑期多校训练营2——Product of GCDs

Product of GCDs

题意：

求给定的多重集合的 求所有长度为k的子集 的$gcd$的 乘积。

思路：

m x = m a x { m u l t i p l e _ s e t k } a n s = ∏ i = 1 m x i k i \begin{aligned} mx&=max\{multiple\_set_k\}\\ ans&=\prod_{i=1}^{mx}i^{k_i}\\ \end{aligned}\\ mxans​=max{multiple_setk​}=i=1∏mx​iki​​
$k_i$表示在集合中的长度为$k$且$gcd$为$i$的子集的个数。
先求出$i$的倍数个数$num$，再求出$dp[i]=\left(\genfrac{}{}{0ex}{}{num}{k}\right)$，相当于在$i$的倍数中取$k$个数，因为其中包含$gcd$为$2i，3i，4i...$的情况，所以要减到这些集合产生的贡献$dp[i]=dp[i]-dp[2i]-dp[3i]...$。(容斥一下)

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N = 1e7+10;
ll n, m, mod;
ll read() {
    ll x = 0; ll f = 1; char s = getchar();
    while(s < '0' || s > '9') {if(s == '-') f = -1; s = getchar();}
    while(s >= '0' && s <= '9') {x = (x << 3) + (x << 1) + s - 48; s = getchar();}
    return x * f;
}
inline ll multi(ll x, ll y, ll mod) {
    ll tmp=(x*y-(ll)((long double)x/mod*y+1.0e-8)*mod);
    return tmp<0 ? tmp+mod : tmp;
}
int p[N], tot;
bool st[N];
void init() {   
    for(int i=2; i<N; i++) {
        if(!st[i]) p[tot++] = i;
        for(int j=0; j<tot&&1ll*i*p[j]<N; j++) {
            st[1ll*i*p[j]] = 1;
            if(i % p[j] == 0) break;
        }
    }
}

ll qpow(ll x, ll y, ll mod) {
    ll ans = 1;
    while(y) {
        if(y & 1) ans = multi(ans, x, mod);
        x = multi(x, x, mod);
        y >>= 1;
    }
    return ans;
}
const int M = 8e5+10;
ll c[M][35];
void Init(ll n, ll mod) {
    c[0][0] = 1;
    for(int i=1; i<=n; i++) {
        for(int j=0; j<=min(31, i); j++) {
            if (!j) {
                c[i][j] = 1;
            }
            else {
                c[i][j] = c[i-1][j] + c[i-1][j-1];
                while(c[i][j]-mod >= mod) c[i][j] -= mod;
            }
        }
    }
}
ll C(ll n, ll m, ll mod) {
    if(n < m) return 0;
    if(m == 0 || n == 0) return 1;
    return c[n][m];
}
ll s[N];

bool Miller_Rabin(ll x) {
	if(x == 2 || x == 3) return true; 
	if(!(x&1)) return false;
	
	ll s = 0, t = x-1;
	while(!(t & 1)) s++, t >>= 1;
	
	for(int i=0; i<5&&p[i]<x; i++) {
		ll a = p[i], k;
		ll b = qpow(a, t, x);// 计算a^t
		for(int j=0; j<s; j++) {
			k = multi(b, b, x); //计算 a^t * a^t, .....
			if(k == 1 && b != 1 && b != x-1) return false;//二次探测（注意二次探测成立得条件之一是方程等于1 <=> k==1）
			b = k;
		}
		if(b != 1) return false; //验证费马定理
	}
	return true;
}
ll cnt[N], dp[N];
ll Phi(ll x) {
    ll ans = x;
    for(int i=0; i<tot&&1ll*p[i]*p[i]<=x; i++) {
        if(x % p[i] == 0) {
            ans = ans / p[i] * (p[i] - 1);
            while(x % p[i] == 0) x /= p[i];
        }
    }
    if(x > 1) ans = ans / x * (x - 1);
    return ans;
}
int main() {
#ifndef ONLINE_JUDGE
    freopen("in.txt", "r", stdin);
    freopen("out.txt", "w", stdout);
    int StartTime = clock();
#endif
	init();
	ll t;
	t = read();
	while(t--) {
		n = read(), m = read(), mod = read();
		ll phi, x, mx = 0, ans = 1;
		if(Miller_Rabin(mod)) phi = mod-1;//mod可能很大，用Miller_Rabin判断一下
		else phi = Phi(mod);//套定义求phi
		Init(n, phi);
		for(int i=1; i<=n; i++) {
			x = read();
			cnt[x]++;
			mx = max(x, mx);
		}
		for(int i = mx; i>=1; i--) {
			//求k_i
			ll now = 0;
			for(int j=i; j<=mx; j+=i) {
				now += cnt[j];
			}
			dp[i] = C(now, m, phi);
			for(int j=i+i; j<=mx; j+=i) {
				dp[i] -= dp[j];
			}
			
			//扩展欧拉定理
			while(dp[i]-phi >= phi) dp[i] -= phi;
			if(dp[i] < 0) dp[i] = (dp[i] % phi + phi) % phi;
			
			ans = multi(ans, qpow(i, dp[i], mod), mod);
		}
		printf("%lld\n", ans%mod);
		for(int i=0; i<=mx; i++) cnt[i] = dp[i] = 0;
	}

#ifndef ONLINE_JUDGE
    printf("Run_Time = %d ms\n", clock() - StartTime);
#endif

    return 0;
}