HDU 4675 GCD of Sequence（莫比乌斯反演）

http://acm.hdu.edu.cn/showproblem.php?pid=4675
题意：给 N $N$ 个数 ai $a_i$，让 ai $a_i$ 变成 bi $b_i$，变的范围为 1 $1$~ M $M$，但是要求要有 K $K$ 个位置的 ai和bi $a_i和b_i$ 不相等，要求输出 b $b$ 这个数组的所有元素的 gcd=1,2,3....M $gcd=1,2,\mathrm{3....}M$ 的方案数

假设我们要求 gcd=d $gcd=d$ 的倍数的方案数：
现在看这 N $N$ 个 ai $a_i$ 中，有 c1 $c1$ 个数是 d $d$ 的倍数，剩下的 c2 $c2$ 个不是 d $d$ 的倍数（ c1+c2=N $c1+c2=N$）

看一个例子 a=2,4,6,5,7 $a=2,4,6,5,7$，现在求 d=2 $d=2$ 的倍数的个数， c1=3,c2=2,K=3,M=10 $c1=3,c2=2,K=3,M=10$

这里因为 5和7 $5和7$ 不是 2 $2$ 的倍数，所以他们肯定要变，而且变了之后就与原来的对应的 ai $a_i$ 不相等了，所以在 K=3 $K=3$ 个不相等的数中就占了两个，而且每个数都能换成 { 2,4,6,8,10 $2,4,6,8,10$} 这 5 $5$ 个数，其实就是 Md $\frac{M}{d}$ 这么多个因此这里的 ans1=(Md)c2 $ans1=(\frac{M}{d}{)}^{c_2}$ 个

然后就是这 c1 $c_1$ 个里面还要选一个换，就是 C1c1 $C_{c1}^1$。但是这里和上面就不一样了，不能换到和自己一样的，就只有 Md−1 $\frac{M}{d}-1$ 个，因此这里的 ans2=CK−c2c1∗(Md−1)K−c2 $ans2=C_{c1}^{K-c_2}\ast (\frac{M}{d}-1{)}^{K-c2}$

其实这个就是求的 F(d) $F\left(d\right)$，然后反演计算 f(d) $f\left(d\right)$ 就行了

一下代码会超时，我试了一下，不把 f(d) $f\left(d\right)$ 以及 F(d) $F\left(d\right)$ 写成函数就能过，但是写成函数思路很清晰

#include"bits/stdc++.h"
#define C(n,m) ((long long)fac[(n)]*inv[(m)]%MOD*inv[(n)-(m)]%MOD)
using namespace std;
typedef long long LL;
const int maxn=3e5+5;
const int MOD=1e9+7;
LL a[maxn],cnt[maxn];//cnt[i]表示是i的倍数的个数 
char mu[maxn];
vector<int>prime;
bool vis[maxn];
LL fac[maxn]={1,1},inv[maxn]={1,1};
int N,M,K;
LL Ans[maxn];
LL Mp[maxn];//记忆化 
long long ksm(long long a,long long b,long long mod)
{
    long long res=1,base=a;
    while(b)
    {
        if(b&1)res=(res*base)%mod;
        base=(base*base)%mod;
        b>>=1;
    }
    return res;
}

void Init(int NN)
{
    memset(vis,1,sizeof(vis));
    mu[1]=1;
    for(int i=2;i<=NN;i++)
    {

        if(vis[i])
        {
            prime.push_back(i);
            mu[i]=-1;
        }
        for(int j=0;j<prime.size()&&i*prime[j]<=NN;j++)
        {
            vis[i*prime[j]]=0;
            if(i%prime[j]==0)
            {
                mu[i*prime[j]]=0;
                break;
            }
            else mu[i*prime[j]]=-mu[i];
        }
        fac[i]=(long long)fac[i-1]*i%MOD;
        inv[i]=ksm(fac[i],MOD-2,MOD);
    }
}
long long F(int d)//F(d)
{
    if(Mp[d])return Mp[d];
    LL c1=cnt[d],c2=N-cnt[d];
    long long ans1=ksm(M/d,c2,MOD);
    long long ans2=(long long)C(c1,K-c2)*ksm(M/d-1,K-c2,MOD)%MOD;
    Mp[d]=ans1*ans2%MOD;
    return Mp[d];
}
long long f(int x)//f(d)
{
    if(N-cnt[x]>K)return 0LL;
    long long res=0;
    for(int d=x;d<=M;d+=x)
    {
        res+=mu[d/x]*F(d);
        res=(res+MOD)%MOD;
    }
    return res;
}

int main()
{
    Init(300000);

    while(cin>>N>>M>>K)
    {
        memset(cnt,0,sizeof(cnt));
        for(int i=1;i<=N;i++)
        {
            scanf("%lld",&a[i]); 
            cnt[a[i]]++;
        }
        for(int i=1;i<=M;i++)
        {
            for(int j=i*2;j<=M;j+=i)
            {
                cnt[i]+=cnt[j];             //求每个数的倍数的个数 
            }
        }
        for(int i=M;i>=1;i--)Ans[i]=f(i);   //之所以倒起来是因为能够用到之前的信息，f(1)计算量最大 
        for(int i=1;i<=M;i++)printf("%lld ",Ans[i]);
        cout<<endl;
        memset(Mp,0,sizeof(Mp));
    }
}

以下是不调用函数，就能过。。。

#include"bits/stdc++.h"
#define C(n,m) ((long long)fac[(n)]*inv[(m)]%MOD*inv[(n)-(m)]%MOD)
using namespace std;
typedef long long LL;
const int maxn=3e5+5;
const int MOD=1e9+7;
LL a[maxn],cnt[maxn];//cnt[i]表示是i的倍数的个数 
char mu[maxn];
vector<int>prime;
bool vis[maxn];
LL fac[maxn]={1,1},inv[maxn]={1,1};
int N,M,K;
LL Ans[maxn];
LL Mp[maxn];//记忆化 
long long ksm(long long a,long long b,long long mod)
{
    long long res=1,base=a;
    while(b)
    {
        if(b&1)res=(res*base)%mod;
        base=(base*base)%mod;
        b>>=1;
    }
    return res;
}

void Init(int NN)
{
    memset(vis,1,sizeof(vis));
    mu[1]=1;
    for(int i=2;i<=NN;i++)
    {

        if(vis[i])
        {
            prime.push_back(i);
            mu[i]=-1;
        }
        for(int j=0;j<prime.size()&&i*prime[j]<=NN;j++)
        {
            vis[i*prime[j]]=0;
            if(i%prime[j]==0)
            {
                mu[i*prime[j]]=0;
                break;
            }
            else mu[i*prime[j]]=-mu[i];
        }
        fac[i]=(long long)fac[i-1]*i%MOD;
        inv[i]=ksm(fac[i],MOD-2,MOD);
    }
}
long long F(int d)//F(d)
{
    if(Mp[d])return Mp[d];
    LL c1=cnt[d],c2=N-cnt[d];
    long long ans1=ksm(M/d,c2,MOD);

    long long ans2=(long long)C(c1,K-c2)*ksm(M/d-1,K-c2,MOD)%MOD;
    Mp[d]=ans1*ans2%MOD;
    return Mp[d];
}
long long f(int x)//f(d)
{
    if(N-cnt[x]>K)return 0LL;
    long long res=0;
    for(int d=x;d<=M;d+=x)
    {
        res+=(LL)mu[d/x]*F(d);
        res=(res+MOD)%MOD;
    }
    return res;
}

int main()
{
    Init(300000);

    while(cin>>N>>M>>K)
    {
        memset(cnt,0,sizeof(cnt));
        for(int i=1;i<=N;i++)
        {
            scanf("%lld",&a[i]); 
            cnt[a[i]]++;
        }
        for(int i=1;i<=M;i++)
        {
            for(int j=i*2;j<=M;j+=i)
            {
                cnt[i]+=cnt[j];             //求每个数的倍数的个数 
            }
        }
        for(int i=M;i>=1;i--)
        {
            if(N-cnt[i]>K)
            {
                Ans[i]=0;
                continue;
            }
            /*这一段相当于F(d)函数,Mp[d]就是F[d]*/
            LL c1=cnt[i],c2=N-cnt[i];
            long long ans1=ksm(M/i,c2,MOD);
            long long ans2=(long long)C(c1,K-c2)*ksm(M/i-1,K-c2,MOD)%MOD;
            Mp[i]=ans1*ans2%MOD;
            /*这一段相当于F(d)函数*/

            /*这一段相当于f(d)函数*/
            long long res=0;
            for(int d=i;d<=M;d+=i)
            {
                res+=(LL)mu[d/i]*Mp[d];
                res=(res+MOD)%MOD;
            }
            /*这一段相当于f(d)函数*/


            Ans[i]=res; 
        }
        for(int i=1;i<=M;i++)
        {
            if(i<M)printf("%lld ",Ans[i]);
            else printf("%lld\n",Ans[i]);
        }
        memset(Mp,0,sizeof(Mp));
    }
}

然后就是别人直接容斥~莫比乌斯函数都懒得弄了，很厉害~