题目：https://ac.nowcoder.com/acm/contest/11253/J
大意为给一个数组a，求数组a中所有大小为k的子集的gcd的乘积。
根据惯例，估计是转化为枚举gcd后统计贡献。对于每个i，我们需要求有几个子集的gcd为i，考虑组合数为，num为数组中i的倍数的个数，但这样统计会包含gcd为i的倍数的集合，因此，考虑从大到小枚举，枚举一个i后，减去gcd为i的倍数的集合的个数，即可得到gcd为i的集合的个数。
综上所述，假设gcd为i的集合的个数为，则答案为
。
但是，还没完，因为数字有这么大，可想而知
的值一定不小，因此，考虑欧拉降幂。我们可以先通过扩展欧拉定理求出
（或者是
）。由于
为一堆组合数的累加，因此，考虑递推预处理组合数（不能用逆元，因为p不一定是质数）。
那么，对于高达1e14的p求，自然是要分解质因数后容斥的。
所以良心出题人为什么不直接开1e18，这样大家就都用pr了
题解说常熟小的暴力可以卡过去，但我是不行，因此，分解质因数用pr解决。
P.S. 好像很多同学关心加不加phi的问题，，我发现欧拉降幂挺多题其实都没故意卡这点。。偷懒点就不写了吧
#include <bits/stdc++.h>
#define inf 0x3f3f3f3f
#define IOS ios_base::sync_with_stdio(0); cin.tie(0);
#define rep(i, a, n) for(int i = a; i <= n; ++ i)
#define per(i, a, n) for(int i = n; i >= a; -- i)
//#define ONLINE_JUDGE
using namespace std;
typedef long long ll;
const int mod=1e9+7;
template<typename T>void write(T x)
{
if(x<0)
{
putchar('-');
x=-x;
}
if(x>9)
{
write(x/10);
}
putchar(x%10+'0');
}
namespace FastIO {
const int SIZE = 1 << 16;
char buf[SIZE], str[64];
int l = SIZE, r = SIZE;
int read(char *s) {
while (r) {
for (; l < r && buf[l] <= ' '; l++);
if (l < r) break;
l = 0, r = int(fread(buf, 1, SIZE, stdin));
}
int cur = 0;
while (r) {
for (; l < r && buf[l] > ' '; l++) s[cur++] = buf[l];
if (l < r) break;
l = 0, r = int(fread(buf, 1, SIZE, stdin));
}
s[cur] = '\0';
return cur;
}
template<typename type>
bool read(type &x, int len = 0, int cur = 0, bool flag = false) {
if (!(len = read(str))) return false;
if (str[cur] == '-') flag = true, cur++;
for (x = 0; cur < len; cur++) x = x * 10 + str[cur] - '0';
if (flag) x = -x;
return true;
}
template <typename type>
type read(int len = 0, int cur = 0, bool flag = false, type x = 0) {
if (!(len = read(str))) return false;
if (str[cur] == '-') flag = true, cur++;
for (x = 0; cur < len; cur++) x = x * 10 + str[cur] - '0';
return flag ? -x : x;
}
} using FastIO::read;
const int maxn=4e4+10;
const int maxai=8e4+10;
const int maxlen=1e7+10;
int t,n,k;
ll p;
int a[maxn];
int cnt[maxai];
ll dp[maxai];
inline ll mul(ll a,ll b,ll mod){
return ((a * b - (long long)((long long)((long double)a / mod * b + 1e-3) * mod)) % mod + mod) % mod;
}
inline ll ksm(ll a,ll n,ll mod){
ll ans=1;
assert(n>=0);
while(n){
if(n&1) ans=mul(ans,a,mod);
a=mul(a,a,mod);
n>>=1;
}
return ans%mod;
}
namespace PHR{
typedef long long LL;
ll n, Max = 0;
inline ll quickpow(ll a, ll b, ll mod)
{
ll ret = 1;
for(; b; b >>= 1, a = mul(a, a, mod))
if(b & 1) ret = mul(ret, a, mod);
return ret;
}
inline bool test(ll a, ll d, ll n)
{
ll t = quickpow(a, d, n);
if(t == 1) return 1;
while(d != n - 1 && t != n - 1 && t != 1) t = mul(t, t, n), d <<= 1;
return t == n - 1; //这里就不用判1了，因为只可能在while前出现1的情况
}
int a[] = {2, 3, 5, 7, 11};
inline bool miller_rabin(ll n)
{
if(n == 1) return 0;
for(int i = 0; i < 5; ++i) //先粗筛一遍
{
if(n == a[i]) return 1;
if(!(n % a[i])) return 0;
}
ll d = n - 1;
while(!(d & 1)) d >>= 1;
for(int i = 1; i <= 5; ++i) //搞五遍
{
ll a = rand() % (n - 2) + 2;//x属于[2, n - 1]
if(!test(a, d, n)) return 0;
}
return 1;
}
inline ll gcd(ll a, ll b) {return b ? gcd(b, a % b) : a;}
inline ll f(ll x, ll a, ll mod) {return (mul(x, x, mod) + a) % mod;}
const int M = (1 << 7) - 1; //我也不知道为啥M是这个数……
ll pollard_rho(ll n) //倍增版！减少gcd调用次数。（好像不用判环？）
{
for(int i = 0; i < 5; ++i) if(n % a[i] == 0) return a[i];
ll x = rand(), y = x, t = 1, a = rand() % (n - 2) + 2;
for(int k = 2;; k <<= 1, y = x)
{
ll q = 1;
for(int i = 1; i <= k; ++i)
{
x = f(x, a, n);
q = mul(q, abs(x - y), n);
// if(i >= M) //不等价！
if(!(i & M)) //超过了2 ^ 7，再用gcd
{
t = gcd(q, n);
if(t > 1) break; //找到了
}
}
if(t > 1 || (t = gcd(q, n)) > 1) break; //之所以再求一遍，是处理刚开始k < M的情况
}
return t;
}
LL factor[100],total=0;
inline void find(ll x)
{
if(x == 1 || x <= Max) return;
if(miller_rabin(x)) {factor[++total]=x; return;}
ll p = x;
while(p == x) p = pollard_rho(x);
while(x % p == 0) x /= p;
find(p); find(x);
}
}
ll x;
ll c[maxai][32];
ll PHIP;
inline ll C(ll n,ll k){
if(n<k)return 0;
return c[n][k];
}
inline ll phi(ll x){
ll ans=x;
for(ll i=2;1ll*i*i<=x;++i){
if(x%i==0){
ans=ans/i*(i-1);
while(x%i==0)x/=i;
}
}
if(x>1)ans=ans/x*(x-1);
return ans;
}
signed main()
{
#ifndef ONLINE_JUDGE
freopen("in.txt","r",stdin);
freopen("out.txt","w",stdout);
#endif
//clock_t c1 = clock();
//===========================================================
read(t);
//init();
while(t--){
read(n),read(k),read(p);
PHR::total=0;
PHR::find(p);
PHIP=p;
for(int i=1;i<=PHR::total;++i)PHIP=PHIP/PHR::factor[i]*(PHR::factor[i]-1);
c[0][0]=1;
int mx=0;
for(int i=1;i<=n;++i){
int x;read(x);
++cnt[x];
mx=max(mx,x);
}
for(int i=1;i<=n;++i){
c[i][0]=1;;
if(i<32)c[i][i]=1;
for(int j=1;j<i&&j<=31;++j){
c[i][j]=(c[i-1][j]+c[i-1][j-1]);
if(c[i][j]>=PHIP)c[i][j]-=PHIP;
}
}
ll ans=1;
for(int i=mx;i>=1;--i){
int now=0;
for(int j=i;j<=mx;j+=i){
now+=cnt[j];
}
dp[i]=C(now,k);
for(int j=i+i;j<=mx;j+=i){
dp[i]-=dp[j];
}
dp[i]=(dp[i]%PHIP+PHIP)%PHIP;
ans=mul(ans,ksm(i,dp[i],p),p);
}
for(int i=0;i<=mx;++i){
cnt[i]=dp[i]=0;
}
printf("%lld\n",ans);
}
//===========================================================
//std::cerr << "Time:" << clock() - c1 << "ms" << std::endl;
return 0;
}
水群的时候看到有聚聚给了个hack样例，扩展欧拉果然挂了，这里再补个修改后的代码：
#include <bits/stdc++.h>
#define inf 0x3f3f3f3f
#define IOS ios_base::sync_with_stdio(0); cin.tie(0);
#define rep(i, a, n) for(int i = a; i <= n; ++ i)
#define per(i, a, n) for(int i = n; i >= a; -- i)
//#define ONLINE_JUDGE
using namespace std;
typedef long long ll;
const int mod=1e9+7;
template<typename T>void write(T x)
{
if(x<0)
{
putchar('-');
x=-x;
}
if(x>9)
{
write(x/10);
}
putchar(x%10+'0');
}
namespace FastIO {
const int SIZE = 1 << 16;
char buf[SIZE], str[64];
int l = SIZE, r = SIZE;
int read(char *s) {
while (r) {
for (; l < r && buf[l] <= ' '; l++);
if (l < r) break;
l = 0, r = int(fread(buf, 1, SIZE, stdin));
}
int cur = 0;
while (r) {
for (; l < r && buf[l] > ' '; l++) s[cur++] = buf[l];
if (l < r) break;
l = 0, r = int(fread(buf, 1, SIZE, stdin));
}
s[cur] = '\0';
return cur;
}
template<typename type>
bool read(type &x, int len = 0, int cur = 0, bool flag = false) {
if (!(len = read(str))) return false;
if (str[cur] == '-') flag = true, cur++;
for (x = 0; cur < len; cur++) x = x * 10 + str[cur] - '0';
if (flag) x = -x;
return true;
}
template <typename type>
type read(int len = 0, int cur = 0, bool flag = false, type x = 0) {
if (!(len = read(str))) return false;
if (str[cur] == '-') flag = true, cur++;
for (x = 0; cur < len; cur++) x = x * 10 + str[cur] - '0';
return flag ? -x : x;
}
} using FastIO::read;
const int maxn=4e4+10;
const int maxai=8e4+10;
const int maxlen=1e7+10;
int t,n,k;
ll p;
int a[maxn];
int cnt[maxai];
ll dp[maxai];
inline ll mul(ll a,ll b,ll mod){
return ((a * b - (long long)((long long)((long double)a / mod * b + 1e-3) * mod)) % mod + mod) % mod;
}
inline ll ksm(ll a,ll n,ll mod){
ll ans=1;
assert(n>=0);
while(n){
if(n&1) ans=mul(ans,a,mod);
a=mul(a,a,mod);
n>>=1;
}
return ans%mod;
}
namespace PHR{
typedef long long LL;
ll n, Max = 0;
inline ll quickpow(ll a, ll b, ll mod)
{
ll ret = 1;
for(; b; b >>= 1, a = mul(a, a, mod))
if(b & 1) ret = mul(ret, a, mod);
return ret;
}
inline bool test(ll a, ll d, ll n)
{
ll t = quickpow(a, d, n);
if(t == 1) return 1;
while(d != n - 1 && t != n - 1 && t != 1) t = mul(t, t, n), d <<= 1;
return t == n - 1; //这里就不用判1了，因为只可能在while前出现1的情况
}
int a[] = {2, 3, 5, 7, 11};
inline bool miller_rabin(ll n)
{
if(n == 1) return 0;
for(int i = 0; i < 5; ++i) //先粗筛一遍
{
if(n == a[i]) return 1;
if(!(n % a[i])) return 0;
}
ll d = n - 1;
while(!(d & 1)) d >>= 1;
for(int i = 1; i <= 5; ++i) //搞五遍
{
ll a = rand() % (n - 2) + 2;//x属于[2, n - 1]
if(!test(a, d, n)) return 0;
}
return 1;
}
inline ll gcd(ll a, ll b) {return b ? gcd(b, a % b) : a;}
inline ll f(ll x, ll a, ll mod) {return (mul(x, x, mod) + a) % mod;}
const int M = (1 << 7) - 1; //我也不知道为啥M是这个数……
ll pollard_rho(ll n) //倍增版！减少gcd调用次数。（好像不用判环？）
{
for(int i = 0; i < 5; ++i) if(n % a[i] == 0) return a[i];
ll x = rand(), y = x, t = 1, a = rand() % (n - 2) + 2;
for(int k = 2;; k <<= 1, y = x)
{
ll q = 1;
for(int i = 1; i <= k; ++i)
{
x = f(x, a, n);
q = mul(q, abs(x - y), n);
// if(i >= M) //不等价！
if(!(i & M)) //超过了2 ^ 7，再用gcd
{
t = gcd(q, n);
if(t > 1) break; //找到了
}
}
if(t > 1 || (t = gcd(q, n)) > 1) break; //之所以再求一遍，是处理刚开始k < M的情况
}
return t;
}
LL factor[100],total=0;
inline void find(ll x)
{
if(x == 1 || x <= Max) return;
if(miller_rabin(x)) {factor[++total]=x; return;}
ll p = x;
while(p == x) p = pollard_rho(x);
while(x % p == 0) x /= p;
find(p); find(x);
}
}
ll x;
ll c[maxai][32];
ll PHIP;
inline ll C(ll n,ll k){
if(n<k)return 0;
return c[n][k];
}
inline ll phi(ll x){
ll ans=x;
for(ll i=2;1ll*i*i<=x;++i){
if(x%i==0){
ans=ans/i*(i-1);
while(x%i==0)x/=i;
}
}
if(x>1)ans=ans/x*(x-1);
return ans;
}
signed main()
{
#ifndef ONLINE_JUDGE
freopen("in.txt","r",stdin);
freopen("out.txt","w",stdout);
#endif
//clock_t c1 = clock();
//===========================================================
read(t);
//init();
while(t--){
read(n),read(k),read(p);
PHR::total=0;
PHR::find(p);
PHIP=p;
for(int i=1;i<=PHR::total;++i)PHIP=PHIP/PHR::factor[i]*(PHR::factor[i]-1);
c[0][0]=1;
int mx=0;
for(int i=1;i<=n;++i){
int x;read(x);
++cnt[x];
mx=max(mx,x);
}
for(int i=1;i<=n;++i){
c[i][0]=1;;
if(i<32)c[i][i]=1;
for(int j=1;j<i&&j<=31;++j){
c[i][j]=(c[i-1][j]+c[i-1][j-1]);
while(c[i][j]-PHIP>=PHIP)c[i][j]-=PHIP;
}
}
ll ans=1;
for(int i=mx;i>=1;--i){
int now=0;
for(int j=i;j<=mx;j+=i){
now+=cnt[j];
}
dp[i]=C(now,k);
for(int j=i+i;j<=mx;j+=i){
dp[i]-=dp[j];
}
while(dp[i]-PHIP>=PHIP)dp[i]-=PHIP;
//dp[i]=(dp[i]%PHIP+PHIP)%PHIP;
if(dp[i]<0)dp[i]=(dp[i]%PHIP+PHIP)%PHIP;
ans=mul(ans,ksm(i,dp[i],p),p);
}
for(int i=0;i<=mx;++i){
cnt[i]=dp[i]=0;
}
printf("%lld\n",ans);
}
//===========================================================
//std::cerr << "Time:" << clock() - c1 << "ms" << std::endl;
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
京公网安备 11010502036488号