2021牛客多校第二场 J——Product of GCDs

题目大意

给你 $n$ 个数，让你从中挑出 $k$ 个数组成一个集合，设第 $i$ 个集合的最大公约数是 $g_i$，对于所有集合，求 $\prod g_i$

解题思路

我们可以枚举 $g_i$，显然 $g_i$ 是从 1 到最大值，，设 $g_i$ 的倍数有 $m$ 个数，那么我们算 $g_i$ 的贡献是 $C_m^k$ 吗？，其实不是，因为 $g_i$ 的倍数组成的集合中，可能最大公约数是 $g_i$ 的倍数，所以我们要减去这些值才是正确的。所以最大公约数为 $g_i$ 的贡献为 $dp[i]=C_m^k-C_{m2}^k-C_{m3}^k-...-C_{mn}^k$， 其中 $mn$ 为 $g_i$ 的 $n$ 倍的数量。
所以最终答案为 ${\prod }_{i=1}^{maxx}{i}^{dp[i]}$
首先，这题中的 $p$ 不一定是质数，所以我们不能用逆元的方式求组合数，我们需要递推求组合数。
$$C_n^m=C_{n-1}^m+C_{n-1}^{m-1}$$
记得初始化 $C[0][0]=1$
第二就是 $dp[i]$ 作为幂数可能太大了，我们需要用欧拉降幂处理一下。
第三就是在乘的时候我们需要注意一下可能爆$longlong$，所以我们用快速乘或者 $int128$ 转一下。

Code

#include <bits/stdc++.h>
#define ll long long
#define qc ios::sync_with_stdio(false); cin.tie(0);cout.tie(0)
#define fi first
#define se second
#define PII pair<int, int>
#define PLL pair<ll, ll>
#define pb push_back
using namespace std;
const int MAXN = 2e5 + 7;
const int inf = 0x3f3f3f3f;
const ll INF = 0x3f3f3f3f3f3f3f3f;
const ll mod = 1e9 + 7;
inline int read()
{
    int x=0,f=1;char ch=getchar();
    while (!isdigit(ch)){if (ch=='-') f=-1;ch=getchar();}
    while (isdigit(ch)){x=x*10+ch-48;ch=getchar();}
    return x*f;
}
int n, k;
ll p;
ll c[MAXN][32];
ll C(int n, int m){
	if(m > n)
		return 0;
	return c[n][m];
}
int cnt[MAXN];
ll dp[MAXN];
inline ll mul(ll a,ll b,ll mod){
    return ((a * b - (long long)((long long)((long double)a / mod * b + 1e-3) * mod)) % mod + mod) % mod;
}
ll qpow(ll x, ll y){
	ll ret = 1;
	while(y){
		if(y&1) ret = mul(ret, x, p);
		x = mul(x, x, p);
		y >>= 1;
	}
	return ret;
}
 
inline ll ksm(ll a,ll n,ll mod){
    ll ans=1;
    assert(n>=0);
    while(n){
        if(n&1) ans=mul(ans,a,mod);
        a=mul(a,a,mod);
        n>>=1;
    }
    return ans%mod;
}
namespace PHR{
    typedef long long LL;
    ll n, Max = 0;
    inline ll quickpow(ll a, ll b, ll mod)
    {
        ll ret = 1;
        for(; b; b >>= 1, a = mul(a, a, mod))
            if(b & 1) ret = mul(ret, a, mod);
        return ret;
    }
 
    inline bool test(ll a, ll d, ll n)
    {
        ll t = quickpow(a, d, n);
        if(t == 1) return 1;
        while(d != n - 1 && t != n - 1 && t != 1) t = mul(t, t, n), d <<= 1;
        return t == n - 1;                       //这里就不用判1了，因为只可能在while前出现1的情况
    }
    int a[] = {2, 3, 5, 7, 11};
	// 不是素数返回0 
    inline bool miller_rabin(ll n)
    {
        if(n == 1) return 0;
        for(int i = 0; i < 5; ++i)        //先粗筛一遍
        {
            if(n == a[i]) return 1;
            if(!(n % a[i])) return 0;
        }
        ll d = n - 1;
        while(!(d & 1)) d >>= 1;
        for(int i = 1; i <= 5; ++i)        //搞五遍
        {
            ll a = rand() % (n - 2) + 2;//x属于[2, n - 1]
            if(!test(a, d, n)) return 0;
        }
        return 1;
    }
 
    inline ll gcd(ll a, ll b) {return b ? gcd(b, a % b) : a;}
    inline ll f(ll x, ll a, ll mod) {return (mul(x, x, mod) + a) % mod;}
    const int M = (1 << 7) - 1;            //我也不知道为啥M是这个数……
	// 找到一个因子返回
    ll pollard_rho(ll n)                    //倍增版！减少gcd调用次数。（好像不用判环？）
    {
        for(int i = 0; i < 5; ++i) if(n % a[i] == 0) return a[i];
        ll x = rand(), y = x, t = 1, a = rand() % (n - 2) + 2;
        for(int k = 2;; k <<= 1, y = x)
        {
            ll q = 1;
            for(int i = 1; i <= k; ++i)
            {
                x = f(x, a, n);
                q = mul(q, abs(x - y), n);
    //            if(i >= M)    //不等价！
                if(!(i & M))    //超过了2 ^ 7，再用gcd
                {
                    t = gcd(q, n);
                    if(t > 1) break;    //找到了
                }
            }
            if(t > 1 || (t = gcd(q, n)) > 1) break;    //之所以再求一遍，是处理刚开始k < M的情况
        }
        return t;
    }
    ll factor[100],total=0;
    inline void find(ll x)
    {
        if(x == 1 || x <= Max) return;
        if(miller_rabin(x)) {factor[++total]=x; return;}
        ll p = x;
        while(p == x) p = pollard_rho(x);
        while(x % p == 0) x /= p;
        find(p); find(x);
    }
}
 
ll PHIP;
inline ll phi(ll x){
    ll ans=x;
    for(ll i=2;1ll*i*i<=x;++i){
        if(x%i==0){
            ans=ans/i*(i-1);
            while(x%i==0)x/=i;
        }
    }
    if(x>1)ans=ans/x*(x-1);
    return ans;
}
 
void solve(){
	cin >> n >> k >> p;
	int maxx = 0;
	for (int i = 1; i <= n; ++i){
		int x;
		cin >> x;
		cnt[x]++;
		maxx = max(maxx, x);
	}
	PHR::total = 0;
	PHR::find(p);
	PHIP = p;
	for(int i = 1; i <= PHR::total; ++i)
		PHIP = PHIP/PHR::factor[i]*(PHR::factor[i] - 1);
	c[0][0] = 1;

	for (int i = 1; i <= n; ++i){
	    c[i][0] = 1;
		for (int j = 1; j <= 31 && j <= i; ++j){
		    c[i][j] = c[i-1][j] + c[i-1][j-1];
			while(c[i][j] - PHIP >= PHIP) c[i][j] -= PHIP;
		}
	}
	ll ans = 1;
	for(int i = maxx; i >= 1; --i){
		int tmp = 0;
		for(int j = i; j <= maxx; j += i){
			tmp += cnt[j];
		}
		dp[i] = C(tmp, k);
		for(int j = i + i; j <= maxx; j += i){
			dp[i] -= dp[j];
		}
		while(dp[i] - PHIP >= PHIP) dp[i] -= PHIP;
		if(dp[i] < 0)
			dp[i] = (dp[i]%PHIP + PHIP) % PHIP;
		ans = mul(ans, qpow(i, dp[i]), p);
	}
	cout << ans << endl;
	for(int i = 0; i <= maxx; i++){
		cnt[i] = dp[i] = 0;
	}
}

int main()
{
    #ifdef ONLINE_JUDGE
    #else
       freopen("in.txt", "r", stdin);
       freopen("out.txt", "w", stdout);
    #endif

    qc;
    int T;
    cin >> T;
    // T = 1;
    while(T--){

        solve();
    }
    return 0;
}