“21 天好习惯”第一期－17

Product of GCDs
思路：
在集合$S$中任选$k$个数组成一个新集合，求所有新集合$GCD$的乘积，因为$x_i$比较小，很容易想到(总数不容易求，常见的套路就是求每个数的贡献)枚举$GCD=p$
假设$F[x]$表示$GCD==x$的新集合的个数，$f[x]$表示$x|GCD$($GCD$是$x$的倍数)的新集合的个数，$a[x]$表示因子包含$x$的数的个数，$GCD=p$的贡献就是$p^{F[x]}$。
那么$f[x]=C_{a[x]}^k$，且$f[x]={\sum }_{x|i}F[i]$
如果要计算$6^{F[6]}$，那么y由上面的定义可以确定的是$2^{f[2]}\ast 3^{f[3]}$中已经计算过$6^{F[6]}$了，那么可以考虑变成枚举$GCD$是质数。
但是对于质数$p$，$GCD=p$的贡献就是$p^{F[x]}$；对于质数$p^2$，$GCD=p^2$的贡献就是$p^{2\ast F[x]}$。而$p^{f[x]}=p^{F[x]+f[2x]+...+f[x^2]+...}$，所以需要再加上$p^{f[x^2]}$。同理对于质数$p^3$，需要再加上$p^{f[x^3]}$。其它同理。于是$p$的贡献就是$p^{\sum f[p^i]}$。

因为$n=4e4$，可想$f[x]$会很大，所以需要用到欧拉降幂，可用公式法求出$\phi (P)$，然后递推杨辉三角预处理组合数。
这题貌似会卡常，加法的时候用减法代替取余，反之取余我就T了。乘法的时候，因为模$P$比较大，所以可以用$\_\_int128$或龟速乘。

code:

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int maxn = 8e4 + 7,maxm=1e7+30;
int n,k,m;
ll p,f[maxn/2][35],a[maxn];
vector<int>prime;
bool vis[maxm];
void initial() {
    for (int i = 2; i < maxm ; ++i) {
        if (!vis[i]) prime.emplace_back(i);
        for (int j = 0,size=prime.size(); j < size; ++j) {
            if (i * prime[j] >= maxm) break;
            vis[i * prime[j]] = true;
            if (i % prime[j] == 0) break;
        }
    }
}
inline ll Phi(ll p){
    ll x=p;
    for(int i:prime) {
    	if(i>x) break;
    	if(x%i==0){
            while(x%i==0)x/=i;
            p=p/i*(i-1);
        }
	}
    if(x>1)p=p/x*(x-1);
    return p;
}
ll qpow(ll a,ll b) {
	ll res=1;
	while(b) {
		if(b&1) res=(__int128)res*a%p;
		a=(__int128)a*a%p;
		b>>=1;
	}
	return res;
}
inline void solve() {
	cin>>n>>k>>p;
	ll phi=Phi(p); m=0;
	f[0][0]=1;
	for(int i=1;i<=n;++i) {
		f[i][0]=1;
		for(int j=min(i,k);j;--j) {
			f[i][j]=f[i-1][j]+f[i-1][j-1];
			if(f[i][j]>=phi) f[i][j]-=phi;
		}
	}
	for(int i=1,x;i<=n;++i) {
		cin>>x; a[x]+=1; if(x>m) m=x;
	}
	ll ans(1);
	for(auto &i:prime) {
		if(i>m) break;
		for(int k=i+i;k<=m;k+=i) a[i]+=a[k];
		ll tmp=f[a[i]][k];
		for(ll j=i*i;j<=m;j*=i) {
			for(int k=j+j;k<=m;k+=j) a[j]+=a[k];
			if(a[j]<k) break;
			tmp+=f[a[j]][k]; if(tmp>=phi) tmp-=phi;
		}
		if(tmp) ans=(__int128)ans*qpow(i,tmp)%p;
	}
	cout<<ans<<'\n'; 
 	memset(a,0,(m+1)<<3);
}
int main() {
    ios::sync_with_stdio(false);cin.tie(0),cout.tie(0);
    initial();
    int T;
    cin>>T;
    while(T--) solve();
    return 0;
}