bzoj 4305: 数列的GCD 莫比乌斯反演

Description

给出一个长度为N的数列{a[n]}，1<=a[i]<=M(1<=i<=N)。
现在问题是，对于1到M的每个整数d，有多少个不同的数列b[1], b[2], …, b[N]，满足：
(1)1<=b[i]<=M(1<=i<=N)；
(2)gcd(b[1], b[2], …, b[N])=d；
(3)恰好有K个位置i使得a[i]<>bi
注：gcd(x1,x2,…,xn)为x1, x2, …, xn的最大公约数。
输出答案对1,000,000,007取模的值。
Input

第一行包含3个整数，N，M，K。
第二行包含N个整数：a[1], a[2], …, a[N]。
Output

输出M个整数到一行，第i个整数为当d=i时满足条件的不同数列{b[n]}的数目mod 1,000,000,007的值。
Sample Input

3 3 3

3 3 3
Sample Output

7 1 0
HINT

当d=1，{b[n]}可以为：(1, 1, 1), (1, 1, 2), (1, 2, 1), (1, 2, 2), (2, 1, 1), (2, 1, 2), (2, 2, 1)。

当d=2，{b[n]}可以为：(2, 2, 2)。

当d=3，因为{b[n]}必须要有k个数与{a[n]}不同，所以{b[n]}不能为(3, 3, 3)，满足条件的一个都没有。

对于100%的数据，1<=N,M<=300000, 1<=K<=N, 1<=a[i]<=M。

分析：
考虑满足有$k$个位置不同，$gcd$是$d$的倍数的方案数。
我们先预处理出$num[x]$表示$a[i]$中有多少个是$x$的倍数。
我们可以开一个桶记录$a[i]$中的数，然后枚举$x$的倍数求$num[x]$。
设$f[d]$为满足有$k$个位置不同，$gcd$是$d$的倍数的方案数。
其中，
$$f[d]=\left(\genfrac{}{}{0ex}{}{num[d]}{n-k}\right)\ast (\frac{m}{d}-1{)}^{num[d]-n+k}\ast (\frac{m}{d}{)}^{n-num[d]}$$
从$m$个数中选出$d$的倍数有$\frac{m}{d}$种方案。
要从$num$个中选出$n-k$个相同的数，剩下的$num[i]-n+k$个数必须不同。选的方法考虑组合数，剩下的必须不同说明选法少了一种，其余情况都是$\frac{m}{d}$种。
$ans[i]={\sum }_{i|j}f[j]\ast \mu (j/i)$

代码：

/**************************************************************
    Problem: 4305
    User: ypxrain
    Language: C++
    Result: Accepted
    Time:812 ms
    Memory:13012 kb
****************************************************************/
 
#include <iostream>
#include <cstdio>
#include <cmath>
#define LL long long
 
const int maxn=3e5+7;
const LL mod=1e9+7;
 
using namespace std;
 
int n,m,k,x,cnt;
int num[maxn],a[maxn];
int prime[maxn],not_prime[maxn];
LL jc[maxn],f[maxn],mu[maxn],ans;
 
void getmu(int n)
{
    mu[1]=1;
    for (int i=2;i<=n;i++)
    {
        if (!not_prime[i])
        {
            prime[++cnt]=i;
            mu[i]=mod-1;
        }
        for (int j=1;j<=cnt;j++)
        {
            if (prime[j]*i>n) break;
            not_prime[i*prime[j]]=1;
            if (i%prime[j]==0)
            {
                mu[i*prime[j]]=0;
                break;
            }
            else mu[i*prime[j]]=mod-mu[i];
        }
    }
}
 
LL power(LL x,LL y)
{
    if (y==0) return 1;
    if (y==1) return x;
    LL c=power(x,y/2);
    c=(c*c)%mod;
    if (y%2) c=(c*x)%mod;
    return c;
}
 
LL C(LL n,LL m)
{
    if (n<m) return 0;
    return jc[n]*power(jc[m]*jc[n-m]%mod,mod-2)%mod;
}
 
int main()
{
    scanf("%d%d%d",&n,&m,&k);
    for (int i=1;i<=n;i++)
    {
        scanf("%d",&x);
        a[x]++;
    }   
    for (int i=1;i<=m;i++)
    {
        for (int j=i;j<=m;j+=i) num[i]+=a[j];
    }
    getmu(m);
    jc[0]=1;
    for (int i=1;i<=n;i++) jc[i]=(jc[i-1]*(LL)i)%mod;
    for (int i=1;i<=m;i++)
    {
        if (num[i]>=n-k) f[i]=C(num[i],n-k)*power(m/i-1,num[i]-n+k)%mod*power(m/i,n-num[i])%mod;
    }   
    for (int i=1;i<=m;i++)
    {
        ans=0;
        for (int j=1;j<=m;j++)
        {
            if (i*j>m) break;
            ans=(ans+mu[j]*f[i*j]%mod)%mod;
        }
        printf("%d ",ans);
    }
}