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]}*3^{f[3]}$ 中已经计算过 $6^{F[6]}$ 了，那么可以考虑变成枚举 ${GCD}$ 是质数。
但是对于质数 ${p}$ ， ${GCD=p}$ 的贡献就是 ${p^{F[x]}}$ ；对于质数 ${p^2}$ ， ${GCD=p^2}$ 的贡献就是 ${p^{2*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]}$ 会很大，所以需要用到欧拉降幂，可用公式法求出 ${φ(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;
}