2021牛客暑期多校2 J

题目大意

给出多重数集 ，求其所有大小为 的子集的 之积。共 组数据，要求答案 输出。

不保证 中所有数各不相同，不保证 为质数。。

事后交题时感觉 。

思路

枚举所有可能出现的 ，设 的长度为 的子集个数为 ，那么最终答案 {%raw%}{%endraw%}。我们的任务是求出每一个 。

枚举 时，统计 中 的倍数出现的个数 ，从这 个中数字中选 个组成备选集合 ，共 中方案， 中所有数的 必然是 的倍数，从 个方案中去掉 的方案，剩下的方案数即为 。只要倒叙枚举 ，。

题目要求答案对 取模，考虑到 很大，我们要进行欧拉降幂。需要注意的是，当 ，也可直接用下面式子的第三条。

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

最终要枚举 计算答案，对当前于 ，要枚举 计算 ，总的时间复杂度是 {%raw%}。

代码

#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;
}