HDU4675 GCD of Sequence

设g(n)为最大公约数等于n的方案数（即本题答案），f(d)为公约数含有d的方案数。

则有关系：f(d) = sigma(d|n, g(n))

本题关键即在于如何求f(d)。

如果没有（恰好K个数与序列a中不同）这个条件的话，那么函数f(d)就很好求，f(d) = [M/d] ^ N

但是本题有这个条件（废话），所以我们需要对函数f的求法做一点改进。

先用O(n sqrt n)的时间预处理出序列a中约数含有x的统计和cnt[x]，

假设我们现在要构造一个最大公约数为d的方案，我们就可以知道a中能被d整除的数有cnt[d]个，要使得整个序列的最大公约数为d的话，就必须从这cnt[d]个数中选n-k个作为不动项，在把剩下来的数作为不同项，且也必须能被d整除。

（1）如果cnt[d] < n-k， 那么就找不出n-k个不动项，所以此时没有合理方案。

（2）如果cnt[d] >= n-k，那么就从中找出n-k个作为不动项，有C(cnt[d], n-k)种方法，再把另外cnt[d]-(n-k)个约数含有d的数变为其他d的倍数，有([M/d] -1) ^ (cnt[d]-(n-k))种方法，最后再把剩下的n-cnt[d]个原本不是d的倍数的数变成d的倍数，有[M/d] ^ (n-cnt[d])种方法。综上，根据乘法原理，共有C(cnt[d], n-k) * ([M/d] -1) ^ (cnt[d]-(n-k)) * n-cnt[d]种方案数。

既然求出了f(d)，剩下的就是容斥原理，我们可以用筛法求所有g(n)，算法复杂度为O(n log n)。

也可以用mobius反演得到：g(n) = sigma(n|d, mu(d/n)*f(d))，只算一个g(n)的算法复杂度为O(n)，计算所有g(n)的算法复杂度依然为O(n log n)。

普通筛法：

#include <cstdio>
#include <iostream>
#include <cstring>
#include <cmath>
#include <ctime>
#include <vector>
#include <set>
#include <queue>
#include <stack>
#include <map>
#include <algorithm>
using namespace std;
typedef long long ll;
typedef unsigned long long llu;

const int INF = 0x3f3f3f3f;
const int mod = 1000000007;
const int maxn = 300005;
const double eps = 1e-6;
const double PI = acos(-1.0);

int n, m, k;
ll fac[maxn];
ll ans[maxn];
int cnt[maxn];

ll pow_mod(ll x, int y)
{
    if (x==1) return 1;
    ll ret=1;
    while (y)
    {
        if (y&1) ret=(ret*x)%mod;
        x=(x*x)%mod;
        y>>=1;
    }
    return ret;
}

ll C(int x,int y)
{
    if (y>x) return 0;
    ll ret = ((fac[x]*pow_mod(fac[y],mod-2))%mod)*pow_mod(fac[x-y],mod-2)%mod;
    return ret;
}

int main()
{
    fac[0]=1;
    for (int i=1;i<=300000;i++) fac[i]=(fac[i-1]*i)%mod;
    while (scanf("%d%d%d",&n,&m,&k)==3)
    {
        memset(cnt,0,sizeof(cnt));
        for (int i=0;i<n;i++)
        {
            int a, j;
            scanf("%d",&a);
            for (j=1;j*j<a;j++)
            {
                if (a%j==0)
                {
                    cnt[j]++;
                    cnt[a/j]++;
                }
            }
            if (j*j==a) cnt[j]++;
        }

        for (int i=m;i>=1;i--)
        {
            if (cnt[i]>=n-k)
            {
                ans[i] = pow_mod(m/i, n-cnt[i]) * C(cnt[i], n-k) % mod * pow_mod(m/i-1, k-n+cnt[i]) % mod;
                for (int j=i*2;j<=m;j+=i) ans[i] = (ans[i] - ans[j] + mod) % mod;
            }
            else ans[i] = 0;
        }
        for (int i=1;i<m;i++)
            printf("%I64d ",ans[i]);
        printf("%I64d\n",ans[m]);
    }
    return 0;
}

mobius反演：

#include <cstdio>
#include <iostream>
#include <cstring>
#include <cmath>
#include <ctime>
#include <vector>
#include <set>
#include <queue>
#include <stack>
#include <map>
#include <algorithm>
using namespace std;
typedef long long ll;
typedef unsigned long long llu;

const int INF = 0x3f3f3f3f;
const int mod = 1000000007;
const int maxn = 300005;
const double eps = 1e-6;
const double PI = acos(-1.0);

int n, m, k;
ll fac[maxn];
ll ans[maxn], f[maxn];
int cnt[maxn];

bool notprime[maxn];
int prime[maxn],nprime;
int mu[maxn];
void getmu()
{
    memset(notprime,0,sizeof(notprime));
    nprime=0;
    mu[1]=1;
    for (int i=2;i<maxn;i++)
    {
        if (!notprime[i]) {prime[nprime++]=i; mu[i]=-1;}
        for (int j=0;j<nprime&&i*prime[j]<maxn;j++)
        {
            notprime[i*prime[j]]=1;
            if (i%prime[j]) mu[i*prime[j]]=-mu[i];
            else mu[i*prime[j]]=0;
        }
    }
}

ll pow_mod(ll x, int y)
{
    if (x==1) return 1;
    ll ret=1;
    while (y)
    {
        if (y&1) ret=(ret*x)%mod;
        x=(x*x)%mod;
        y>>=1;
    }
    return ret;
}

ll C(int x,int y)
{
    if (y>x) return 0;
    ll ret = ((fac[x]*pow_mod(fac[y],mod-2))%mod)*pow_mod(fac[x-y],mod-2)%mod;
    return ret;
}

int main()
{
    getmu();
    fac[0]=1;
    for (int i=1;i<=300000;i++) fac[i]=(fac[i-1]*i)%mod;
    while (scanf("%d%d%d",&n,&m,&k)==3)
    {
        memset(cnt,0,sizeof(cnt));
        for (int i=0;i<n;i++)
        {
            int a, j;
            scanf("%d",&a);
            for (j=1;j*j<a;j++)
            {
                if (a%j==0)
                {
                    cnt[j]++;
                    cnt[a/j]++;
                }
            }
            if (j*j==a) cnt[j]++;
        }

        for (int i=m;i>=1;i--)
        {
            if (cnt[i]>=n-k) f[i] = C(cnt[i], n-k) * pow_mod(m/i-1, cnt[i]-n+k) % mod * pow_mod(m/i, n-cnt[i]) % mod;
            else f[i] = 0;
        }

        for (int i=1;i<=m;i++)
        {
            ans[i] = 0;
            for (int j=i;j<=m;j+=i)
                ans[i] = (ans[i] + mu[j/i]*f[j] + mod) % mod; // mu可能为-1
        }

        for (int i=1;i<m;i++)
            printf("%I64d ",ans[i]);
        printf("%I64d\n",ans[m]);
    }
    return 0;
}