Product of GCDs 题解(欧拉降幂+贡献)

题目链接

题目大意

给一个数组\(a\)，求数组\(a\)中所有大小为\(k\)的子集的\(\gcd\)的乘积。

题目思路

是一个比较经典的贡献问题

下面说下注意的点

• 求欧拉函数不能直接暴力for 2 3 4 5这样枚举 要直接枚举质数
• 欧拉降幂注意对于大于\(\phi(x)\)的\(y\)，其结果应该是$y; mod;\phi(x)+\phi(x) $ 尽管本题没卡,代码末尾hack数据
• 由于\(\phi(x)\)的值可以达到\(10^{14}\)的级别，所以要用__int128
• 求贡献时候的细节有点多，要考虑如何去重
• 组合数求解的时候数据范围小可以直接递推，由于不一定是质数，所以如果求逆元的话不能快速幂求解

代码

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
//#define __int128 ll
//typedef pair<int,int> pii;
#define fi first
#define se second
#define debug printf("aaaaaaaaaaa\n");
const int maxn=1e7+5,inf=0x3f3f3f3f,mod=1e9+7;
const ll INF=0x3f3f3f3f3f3f3f3f;
int n,k;
ll p;
ll a[40000+5];
ll c[40000+5][30+5];
ll num[80000+5];
ll prime[maxn],cnt;
bool isprime[maxn];
void getprime(int n){
    for(ll i=2;i<=n;i++){//开ll因为后面要计算i*prime[j]
        if(!isprime[i]){
            prime[++cnt]=i;
        }
        for(int j=1;j<=cnt&&i*prime[j]<=n;j++){
            isprime[i*prime[j]]=1;
            if(i%prime[j]==0) break;
        }
    }
}
ll get(ll n){ // 计算欧拉函数
    ll num=n;
    for(int i=1; prime[i]*prime[i]<=n; ++i)
        if(n%prime[i]==0) {
            num=num/prime[i]*(prime[i]-1);
            while(n%prime[i]==0){
                n/=prime[i];
            }
        }
    if(n>1){
        num=num/n*(n-1);
    }
    return num;
}
ll qpow(ll a,ll b,ll mod){
    __int128 ans=1,base=a;
    while(b){
        if(b&1){
            ans=ans*base%mod;
        }
        base=base*base%mod;
        b=b>>1;
    }
    return ans;
}
int main(){
    getprime(10000000);
    int _;scanf("%d",&_);
    while(_--){
        ll ma=0;
        scanf("%d%d%lld",&n,&k,&p);
        for(int i=1;i<=8e4;i++){
            num[i]=0;
        }
        for(int i=1;i<=n;i++){
            scanf("%lld",&a[i]);
            num[a[i]]++;
            ma=max(ma,a[i]);
        }
        ll mod=get(p);
        c[0][0]=1;
        for(int i=1;i<=n;i++){
            c[i][0]=1;
            for(int j=1;j<=min(i,k);j++){
                c[i][j]=c[i-1][j-1]+c[i-1][j];
                if(c[i][j]>=mod){
                    c[i][j]=c[i][j]%mod+mod;
                }
            }
        }
        ll ans=1;
        for(int i=1;prime[i]<=ma;i++){
            for(ll j=prime[i];j<=ma;j*=prime[i]){
                ll temp=0;
                for(ll k=j;k<=ma;k+=j){
                    temp+=num[k];
                }
                if(temp<k) break;
                ans=(__int128)ans*qpow(prime[i],c[temp][k],p)%p;
            }
        }
        printf("%lld\n",ans);
    }
    return 0;
}
//1
//24 7 1038318
//2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
//正确答案是519160