2021牛客暑期多校2 J

题目大意

给出多重数集 $S={x_1,x_2,\cdots,x_n}$ ，求其所有大小为 $k$ 的子集的 $\gcd$ 之积。共 $T$ 组数据，要求答案 $\bmod P$ 输出。

不保证 $S$ 中所有数各不相同，不保证 $P$ 为质数。 $T≤60,10^6≤P≤10^{14},1≤x_i≤8\times 10^4,1≤|S|≤40000,1≤k≤\min(|S|,30)$ 。

事后交题时感觉 $|S|\le 50000$ 。

思路

枚举所有可能出现的 $\gcd$ ，设 $\gcd=g$ 的长度为 $k$ 的子集个数为 $c_g$ ，那么最终答案 {%raw%} $ans=\prod\limits_{g=1}^{\max\limits_{1\le i\le n}\{x_i\}}g^{c_g}$ {%endraw%}。我们的任务是求出每一个 $c_g$ 。

枚举 $g$ 时，统计 $S$ 中 $g$ 的倍数出现的个数 $num_g=\sum\limits_{i=1}^n[g\mid x_i]$ ，从这 $num_g$ 个中数字中选 $k$ 个组成备选集合 $S'$ ，共 $C_{num_g}^k$ 中方案， $S'$ 中所有数的 $\gcd$ 必然是 $g$ 的倍数，从 $C_{num_g}^{k}$ 个方案中去掉 $\gcd=2g,3g\cdots$ 的方案，剩下的方案数即为 $c_g$ 。只要倒叙枚举 $g$ ， $c_g=C_{num_g}^k-c_{2g}-c_{3g}-\cdots$ 。

题目要求答案对 $P$ 取模，考虑到 $c_g$ 很大，我们要进行欧拉降幂。需要注意的是，当 $\gcd(a,m)=1$ ，也可直接用下面式子的第三条。

$扩展欧拉定理\quad a^b \equiv \begin{cases} a^{b\ \bmod\ \varphi(m)}, &\gcd(a,m) = 1, \\ a^b, &\gcd(a,m)\ne 1, b < \varphi(m),\\ a^{b\ \bmod\ \varphi(m) +\varphi(m)}, &\gcd(a,m)\ne 1, b \ge \varphi(m). \end{cases} \pmod m$

$欧拉函数\quad\varphi(m)=m\times \prod\limits_{p_i\in \mathbb{prime},p_i\mid m}\left(1-\frac{1}{p_i}\right)$

考虑到 $P$ 很大，我们用 Pollard-Rho 算法质因数分解得到 $P$ 的所有质因子，再计算 $\varphi(P)$ 。同时由于最后统计答案多次用到组合数，且组合数是作为指数出现，因此可以预处理组合数；因为 $φ (P)$ 不一定为质数，同时 $k$ 不大，所以我们可以用杨辉三角递推组合数。模乘法需要用龟速乘，~~但是有__int128~~。

最终要枚举 $g\in[1,\max\{x_i\}]$ 计算答案，对当前于 $g$ ，要枚举 $c_{2g},c_{3g},\cdots$ 计算 $c_g$ ，总的时间复杂度是 {%raw%} $O\left(\max\limits_{1\le i\le n}\{x_i\}\ln\left(\max\limits_{1\le i\le n}\{x_i\}\right)\right)$ 。

代码

#include<iostream>
#include<ctime>
#include<algorithm>
#include<cstring>
#include<unordered_set>
using namespace std;
typedef long long ll;
typedef __int128 bll;
int _;
unordered_set<ll>primeFactor;
int n, k;
ll P, phiP;
ll combine[80004][32];
int cnt[80004];
ll c[80004];
// 快速幂
ll QuickPow(ll a, ll b, ll MOD)
{
    a %= MOD;
    ll res = 1 % MOD;
    while (b)
    {
        if (b & 1)
        {
            res = (bll)res * (bll)a % MOD;
        }
        b >>= 1;
        a = (bll)a * (bll)a % MOD;
    }
    return res;
}
// 素数测试
bool MillerRabin(ll n, int repeat)
{
    if (n == 2 || n == 3)
    {
        return true;
    }
    else if (n % 2 == 0 || n == 1)
    {
        return false;
    }
    // 将n-1分解成2^s*d
    ll d = n - 1;
    int s = 0;
    while (!(d & 1))
    {
        ++s;
        d >>= 1;
    }
    while (repeat--)
    {
        // 取随机数[2,n-2]
        ll a = rand() % (n - 3) + 2;
        ll x = QuickPow(a, d, n);
        ll y = 0;
        for (int i = 1;i <= s;i++)
        {
            y = (bll)x * (bll)x % n;
            if (y == 1 && x != 1 && x != n - 1)
            {
                return false;
            }
            x = y;
        }
        if (y != 1)
        {
            return false;
        }
    }
    return true;
}
// 返回最大质因子
ll PollardRho(ll n, ll c)
{
    ll x = rand() % (n - 2) + 1;
    ll y = x, i = 1, k = 2;
    while (1)
    {
        ++i;
        x = (bll)x * (bll)x % n + c + n;
        ll d = abs(__gcd(y - x, n));
        if (1 < d && d < n)
        {
            return d;
        }
        else if (y == x)
        {
            return n;
        }
        else if (i == k)
        {
            y = x;
            k <<= 1;
        }
    }
}
// 分解质因数
void FindPrimeFactor(ll n, ll c)
{
    if (n == 1)
    {
        return;
    }
    if (MillerRabin(n, 50))
    {
        primeFactor.emplace(n);
        return;
    }
    ll p = n;
    while (p >= n)
    {
        p = PollardRho(p, c--);
    }
    FindPrimeFactor(p, c);
    FindPrimeFactor(n / p, c);
}
void init()
{
    // 求P的欧拉函数值
    FindPrimeFactor(P, 103);
    phiP = P;
    for (const ll& p : primeFactor)
    {
        phiP = phiP / p * (p - 1);
    }

    // 预处理组合数
    for (int i = 1;i <= n;i++)
    {
        combine[i][0] = combine[i][i] = 1;
        for (int j = 1;j < i && j <= k;j++)
        {
            // 杨辉三角递推
            combine[i][j] = combine[i - 1][j] + combine[i - 1][j - 1];
            // 欧拉降幂
            if (combine[i][j] >= phiP)
            {
                combine[i][j] = combine[i][j] % phiP + phiP;
            }
        }
    }
}
int main()
{
    srand(time(0));
    for (cin >> _;_;_--)
    {
        scanf("%d %d %lld", &n, &k, &P);
        init();
        int maxX = 0;
        for (int i = 1;i <= n;i++)
        {
            int x;
            scanf("%d", &x);
            maxX = max(maxX, x);
            // 预处理cnt数组 记录x出现的次数
            ++cnt[x];
        }
        ll ans = 1;
        for (int g = maxX;g >= 1;g--)
        {
            int num = 0;
            // 得到g的倍数的个数
            for (int x = g;x <= maxX;x += g)
            {
                num += cnt[x];
            }
            c[g] = combine[num][k];
            // 排除公约数为g的倍数的集合
            for (int x = 2 * g;x <= maxX;x += g)
            {
                c[g] -= c[x];
            }
            if (c[g] >= phiP)
            {
                c[g] = c[g] % phiP + phiP;
            }
            else if (c[g] < 0)
            {
                c[g] = (c[g] % phiP + phiP) % phiP;
            }
            ans = (bll)ans * (bll)QuickPow(g, c[g], P) % P;
        }
        printf("%lld\n", ans);
        memset(cnt, 0, sizeof(int) * (maxX + 1));
        memset(c, 0, sizeof(ll) * (maxX + 1));
        primeFactor.clear();
    }
    return 0;
}