分析
我们设
f
(
d
)
f(d)
f(d)为
gcd
=
d
\gcd=d
gcd=d的数量,
F
(
d
)
F(d)
F(d)为
gcd
\gcd
gcd为
d
d
d或者
d
d
d的倍数的数量
我们有反演(假设n=3):
f
(
d
)
=
∑
d
∣
t
μ
(
t
d
)
F
(
t
)
f(d)=\sum_{d \mid t}\mu(\frac{t}{d})F(t)
f(d)=∑d∣tμ(dt)F(t)
所以这道题难在
F
(
t
)
F(t)
F(t)怎么求:
我们设
n
u
m
[
d
]
num[d]
num[d]为
a
a
a数组中,是
d
d
d的倍数的数目,那剩下的
n
−
n
u
m
[
d
]
n-num[d]
n−num[d]是需要改为
d
d
d的倍数。因为
gcd
\gcd
gcd为
d
d
d
我们要保证有正好
k
k
k个不同,
n
−
n
u
m
[
d
]
n-num[d]
n−num[d]是一定要修改的,那可能还有
k
−
(
n
−
n
u
m
[
d
]
)
k-(n-num[d])
k−(n−num[d])个需要修改
这些的贡献为
C
n
u
m
[
d
]
k
+
n
u
m
[
d
]
−
n
C_{num[d]}^{k+num[d]-n}
Cnum[d]k+num[d]−n
那可以修改的数量应该为
m
d
−
1
\frac{m}{d}-1\quad
dm−1(除去本身这个倍数)
所以这部分的总贡献:
C
n
u
m
[
d
]
k
+
n
u
m
[
d
]
−
n
(
m
d
−
1
)
k
+
n
u
m
[
d
]
−
n
C_{num[d]}^{k+num[d]-n}(\frac{m}{d}-1)^{k+num[d]-n}
Cnum[d]k+num[d]−n(dm−1)k+num[d]−n
还有哪些必须修改的贡献,所以最终:
F
(
d
)
=
C
n
u
m
[
d
]
k
+
n
u
m
[
d
]
−
n
(
m
d
−
1
)
k
+
n
u
m
[
d
]
−
n
(
m
d
)
n
−
n
u
m
[
d
]
F(d)=C_{num[d]}^{k+num[d]-n}(\frac{m}{d}-1)^{k+num[d]-n}(\frac{m}{d})^{n-num[d]}
F(d)=Cnum[d]k+num[d]−n(dm−1)k+num[d]−n(dm)n−num[d]
最后我们枚举 d d d,在枚举 d d d的倍数即可 O ( n l o g n ) O(nlogn) O(nlogn)求解
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
typedef unsigned long long ull;
template <typename T>
void out(T x) { cout << x << endl; }
ll fast_pow(ll a, ll b, ll p) {ll c = 1; while(b) { if(b & 1) c = c * a % p; a = a * a % p; b >>= 1;} return c;}
ll exgcd(ll a, ll b, ll &x, ll &y) { if(!b) {x = 1; y = 0; return a; } ll gcd = exgcd(b, a % b, y, x); y-= a / b * x; return gcd; }
const int N = 3e5 + 6;
const int mod = 1e9 + 7;
int prime[N], mu[N], tot;
ll F[N], num[N], fac[N], inv[N];
bool mark[N];
void init()
{
fac[0] = inv[0] = 1;
for(int i = 1; i < N; i ++)
fac[i] = fac[i - 1] * i % mod;
inv[N - 1] = fast_pow(fac[N - 1], mod - 2, mod);
for(int i = N - 2; i >= 1; i --)
inv[i] = inv[i + 1] * (i + 1) % mod;
}
void get_mu()
{
mark[1] = mark[0] = true;
mu[1] = 1;
tot = 0;
for(int i = 2; i < N; i ++)
{
if(!mark[i])
{
prime[tot ++] = i;
mu[i] = -1;
}
for(int j = 0; j < tot && i * prime[j] < N; j ++)
{
mark[i * prime[j]] = true;
if(i % prime[j] == 0)
break;
mu[i * prime[j]] = -mu[i];
}
}
}
ll C(int n, int m)
{
return fac[n] * inv[m] % mod * inv[n - m] % mod;
}
int main()
{
ios::sync_with_stdio(false);
int n, m, k;
get_mu();
init();
while(cin >> n >> m >> k)
{
memset(num, 0, sizeof(num));
for(int i = 1; i <= n; i ++)
{
int x;
cin >> x;
for(int j = 1; j * j <= x; j ++)
{
if(x % j == 0)
{
num[j] ++;
if(j * j != x)
num[x / j] ++;
}
}
}
for(int i = 1; i <= m; i ++)
{
if(k + num[i] - n < 0)
F[i] = 0;
else
F[i] = C(num[i], k - (n - num[i])) * fast_pow(m / i - 1, k + num[i] - n, mod) % mod * fast_pow(m / i, n - num[i], mod) % mod;
}
for(int i = 1; i <= m; i ++)
{
ll ans = 0;
for(int j = i; j <= m; j += i)
{
ans = (ans + 1ll * mu[j / i] * F[j] % mod + mod) % mod;
}
cout << ans << (i != m ? " " : "\n");
}
}
//system("pause");
}