2021牛客多校 Product of GCDs（拓展欧拉定理+欧拉函数）

题目传送门

大意

给一个数组，求数组中所有大小为k的子集的gcd的乘积。
样例：
5 2 20040220
2 4 8 12 14
这个就是从第二行找到k=2个数，（选出来，然后将这些组的gcd乘积求出来）
2-4：2
2-8：2
2-12：2
2-14：2
4-8：4
4-12：4
4-14：2
8-12：4
8-14：2
12-14：2
乘起来，然后对最后答案20040220取模
结果就是8192

思路

其实就是根据素数，来计算它对整体答案的贡献，来找它的倍数的大小，最后计算贡献：
比如现在给你一组数它的最小因子是$p$，那么来讨论它的倍数$2\ast p,3\ast p.......$，最后得出了m组与p相关的倍数，那么这组数对答案的贡献就是用组合数算出来的$C_m^k$因此，对答案的贡献是$p^{C_m^k}$
但是这还不够，我们现在枚举了$p$，但是对于那些有$p^k$(k为正整数)的因子的元素来说，我们对于他们对最后的$gcd$求小了，因此我们对这些$p$的乘方都要进行上述运算。以此类推。

同时由于我们对答案的贡献可能会很大，因此要用欧拉降幂（这一点不用多说），那么我们在递推求组合数的过程中就可以进行欧拉定理的拓展了。

那么这个题就有一个非常容易超时的点就是：我们在求欧拉函数的时候，应该用素数进行枚举，而不是用普通的常数。

#pragma GCC optimize(2)
#include<iostream>
#include<cstring>
#include<cmath>
using namespace std;
typedef long long ll;
const int N=1e7+8;
int primes[N+8], cnt;
bool st[N+8];
ll c[40008][31];
int a[40005];
int ok[80005];
ll euler;
inline int read()
{
    int x=0,f=1;
    char ch=getchar();
    while(ch<'0'||ch>'9')
    {
        if(ch=='-')
            f=-1;
        ch=getchar();
    }
    while(ch>='0' && ch<='9')
        x=x*10+ch-'0',ch=getchar();
    return x*f;
}
inline ll lread()
{
    ll x=0,f=1;
    char ch=getchar();
    while(ch<'0'||ch>'9')
    {
        if(ch=='-')
            f=-1;
        ch=getchar();
    }
    while(ch>='0' && ch<='9')
        x=x*10+ch-'0',ch=getchar();
    return x*f;
}
/*void init_C(int n){
	for (int i = 0; i <= n; i ++ )
    for (int j = 0; j <= i; j ++ ){
    	if(j>=31)break;
    	if (!j) c[i][j] = 1;
        else c[i][j] = (c[i - 1][j] + c[i - 1][j - 1]);
        if(c[i][j]>euler)c[i][j]%=euler,c[i][j]+=euler;
	}
    return;
}*/
void init_primes(){
	for (int i = 2; i <= N; i ++ )
    {
        if (!st[i]) primes[cnt ++ ] = i;
        for (int j = 0; primes[j] <= N / i; j ++ )
        {
            st[primes[j] * i] = true;
            if (i % primes[j] == 0) break;
        }
    }
    return;
}
ll phi(ll x){//求欧拉函数 
    __int128 res = x;
    for (ll i = 0; primes[i]<= x / primes[i]; i ++ )
        if (x % primes[i] == 0)
        {
            res = res / primes[i] * (primes[i] - 1);
            while (x % primes[i] == 0) x /= primes[i];
        }
    if (x > 1) res = res / x * (x - 1);
    return res;
}
ll qmi(ll m, ll k, ll p)
{
    __int128 res = 1 % p, t = m;
    while (k)
    {
        if (k&1) res = res * t % p;
        t = (__int128)t * t % p;
        k >>= 1;
    }
    return res;
}

int main(){
    std::ios::sync_with_stdio(false);
    cin.tie(0); 
    cout.tie(0);
	init_primes();
	int T;
	//cin>>T;
    T=read();
	while(T--){
		int n,k;
		ll modd;
		//cin>>n>>k>>modd;
        n=read();
        k=read();
        modd=lread();
		memset(ok,0,sizeof ok);
        //for(int i=0;i<=80000;i++)ok[i]=0;
		int maxx=0;
		euler=phi(modd);
        //euler=1000000;
        //init_C(n);
        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]>=euler){
                    c[i][j]=c[i][j]%euler+euler;
                }
            }
        }
		for(int i=1;i<=n;i++){
			//cin>>a[i];
            a[i]=read();
			ok[a[i]]++;
			maxx=max(maxx,a[i]);
		}
		//cout<<c[5][2]<<" "<<c[3][2]<<endl;(__int128)
		ll ans=1;
		for(int i=0;primes[i]<=maxx;i++){
			for(ll j=primes[i];j<=maxx;j*=primes[i]){
				int s=0;
				for(ll _p=j;_p<=maxx;_p+=j)s+=ok[_p];
				if(s<k)break;
				ans=(__int128)ans*qmi(primes[i],c[s][k],modd)%modd;

			}
		}
		cout<<ans<<endl;
	} 
}