ACM学习历程—HDU4675 GCD of Sequence（莫比乌斯）

Description

Alice is playing a game with Bob. 
Alice shows N integers a 1, a 2, …, a N, and M, K. She says each integers 1 ≤ a i ≤ M. 
And now Alice wants to ask for each d = 1 to M, how many different sequences b 1, b 2, …, b N. which satisfies : 
1. For each i = 1…N, 1 ≤ b[i] ≤ M 
2. gcd(b 1, b 2, …, b N) = d 
3. There will be exactly K position i that ai != bi (1 ≤ i ≤ n) 

Alice thinks that the answer will be too large. In order not to annoy Bob, she only wants to know the answer modulo 1000000007.Bob can not solve the problem. Now he asks you for HELP! 
Notes: gcd(x 1, x 2, …, x n) is the greatest common divisor of x 1, x 2, …, x n

  

Input

The input contains several test cases, terminated by EOF. 
The first line of each test contains three integers N, M, K. (1 ≤ N, M ≤ 300000, 1 ≤ K ≤ N) 
The second line contains N integers: a 1, a 2, ..., a n (1 ≤ a i ≤ M) which is original sequence. 

  

Output

For each test contains 1 lines : 
The line contains M integer, the i-th integer is the answer shows above when d is the i-th number. 

  

Sample Input

3 3 33 3 33 5 31 2 3 

  

Sample Output

7 1 059 3 0 1 1 

Hint

 In the first test case : when d = 1, {b} can be : (1, 1, 1) (1, 1, 2) (1, 2, 1) (1, 2, 2) (2, 1, 1) (2, 1, 2) (2, 2, 1) when d = 2, {b} can be : (2, 2, 2) And because {b} must have exactly K number(s) different from {a}, so {b} can't be (3, 3, 3), so Answer = 0 

题目关键在于找出什么东西是可以求的。

之前做过莫比乌斯的题目，一般都是d | gcd的比较好求。

那么考虑对于d | gcd的情况F(d)有多少种？

如果当前数列里面有num[d]个数是d的倍数。

首先由于要改掉k个，那么有n-k个数是不变的，如果这些数有不是d的倍数，那么不满足条件了。

即如果n-k > num[d]，F(d) == 0。

于是对于余下的n-k，可以从num[d]个数中选取，即C(num[d], n-k)。

然后对于被改掉了k个数来说，里面肯定有num[d]-(n-k)个数是d的倍数，则有quickPow(m/d-1,num[d]-(n-k))种选法。最后剩余的n-num[i]个数，则有quickPow(m/i, n-num[i])种选法。最后乘起来就是结果了。需要对于0的情况考虑一下。

于是F(d)有了，而F(d) = sum(f(n)) (d|n)。

这一步网上有把F(d)一道一边，f(d)一到另一边计算的。这里采用莫比乌斯反演。

然后C组合数的求法的话，预处理出所有排列数p[]和逆元invp[]，这里的invp是用前一个invp[i-1]乘上i的逆元计算的，i的逆元又是用O(n)算法先处理出来的。总的复杂度是O(n)。

快速幂复杂度是logn。

然后莫比乌斯的复杂度是nlogn。

最后总的复杂度是nlogn的。

代码：

#include <iostream>
#include <cstdio>
#include <cstdlib>
#include <cmath>
#include <cstring>
#include <algorithm>
#include <set>
#include <map>
#include <queue>
#include <string>
#define LL long long
#define MOD 1000000007

using namespace std;

const int maxN = 300010;
int n, m, k;
int num[maxN];
int sum[maxN];
int p[maxN], inv[maxN], invp[maxN];
int u[maxN], prime[maxN];
bool vis[maxN];
int ans[maxN];

//莫比乌斯反演
//F(n)和f(n)是定义在非负整数集合上的两个函数
//并且满足条件F(n) = sum(f(d)) (d|n)
//那么f(n) = sum(u(d)*F(n/d)) (d|n)
//case 1: if d = 1, u(d) = 1
//case 2: if d = p1*p2...pk, (pi均为互异素数), u(d) = (-1)^k
//case 3: else, u(d) = 0;
//性质1: sum(u(d)) (d|n) = 1 (n == 1) or 0 (n > 1)
//性质2: sum(u(d)/d) (d|n) = phi(n)/n
//另一种形式:F(d) = sum(f(n)) (d|n) => f(d) = sum(u(n/d)*F(n)) (d|n)
//线性筛选求莫比乌斯反演函数代码
void mobius()
{
    memset(vis, false,sizeof(vis));
    u[1] = 1;
    int cnt = 0;
    for(int i = 2; i < maxN; i++)
    {
        if(!vis[i])
        {
            prime[cnt++] = i;
            u[i] = -1;
        }
        for(int j = 0; j < cnt && i*prime[j] < maxN; j++)
        {
            vis[i*prime[j]] = true;
            if(i%prime[j])
                u[i*prime[j]] = -u[i];
            else
            {
                u[i*prime[j]] = 0;
                break;
            }
        }
    }
}

void init()
{
    mobius();
    //***预处理所有ｉ在质数MOD下的逆元
    inv[1] = 1;
    for (int i = 2; i < maxN; i++)
        inv[i] = (LL)inv[MOD%i]*(MOD-MOD/i) % MOD;
    p[0] = 1;
    invp[0] = inv[1];
    for (int i = 1; i < maxN; ++i)
    {
        p[i] = (LL)p[i-1]*i%MOD;
        invp[i] = (LL)invp[i-1]*inv[i]%MOD;
    }
}

void input()
{
    memset(num, 0, sizeof(num));
    int tmp;
    for (int i = 0; i < n; ++i)
    {
        scanf("%d", &tmp);
        num[tmp]++;
    }

    for (int i = 1; i <= m; ++i)
    {
        for (int j = i*2; j <= m; j += i)
            num[i] += num[j];
    }
}

inline int C(int x, int y)
{
    int ans;
    ans = (LL)p[x]*invp[x-y]%MOD;
    ans = (LL)ans*invp[y]%MOD;
    return ans;
}

//快速幂m^n
LL quickPow(LL x, int n)
{
    if (n == 0)
        return 1;
    if (x == 0)
        return 0;
    LL a = 1;
    while (n)
    {
        a *= n&1 ? x : 1;
        a %= MOD;
        n >>= 1 ;
        x *= x;
        x %= MOD;
    }
    return a;
}

void work()
{
    memset(sum, 0, sizeof(sum));
    memset(ans, 0, sizeof(ans));
    for (int i = 1; i <= m; ++i)
    {

        if (num[i] < n-k)
            sum[i] = 0;
        else
        {
            sum[i] = C(num[i], n-k)*quickPow(m/i, n-num[i])%MOD;
            sum[i] = quickPow(m/i-1, num[i]-(n-k))*sum[i]%MOD;
        }
    }
    for (int i = 1; i <= m; ++i)
    {
        for (int j = i; j <= m; j += i)
        {
            ans[i] += (LL)u[j/i]*sum[j]%MOD;
            ans[i] = (ans[i]%MOD+MOD)%MOD;
        }
    }
    for (int i = 1; i <= m; ++i)
    {
        if (i != 1)
            printf(" ");
        printf("%d", ans[i]);
    }
    printf("
");
}

int main()
{
    //freopen("test.in", "r", stdin);
    init();
    while (scanf("%d%d%d", &n, &m, &k) != EOF)
    {
        input();
        work();
    }
    return 0;
}