hdu6053-容斥+莫比乌斯函数+优化

题目链接：http://acm.hdu.edu.cn/showproblem.php?pid=6053

题意：给一个序列A，要求构造序列B，使得 Bi <= Ai, gcd(Bi) > 1, 1 <= i <= n, 输出构造的方法数取模1e9+7

思路：

求出序列A的最大值Max，最小值Min

枚举整个序列的{d| 2 <= d <= Min}, 对于每个d来说，每个Bi可以构造的数为d, 2*d, ..., (Ai/d)*d,

所以每个d对最后结果的贡献值为 π(Ai/d), 如果不考虑重复，答案为Σd*π(Ai/d),但是最后会有重复，

例如d = 2时, 3*d = 6, d = 3时, 2*d = 6, d = 6时, 1*d = 6，这样6这个构造方案就重复了，需要减去，

而直接减去d=6时的方案，会导致减多了，需要加上多减的一部分，到这里其实就是容斥了...

对d进行质因数分解，d = a1^p1 * a2^p2 * ... * am^pm，若存在pi>1，则此时的方案数忽略掉，不加也不减，

若m为奇数，即d由奇数个质因数组成，则此时的方案数要加上，若m为偶数，则此时的方案数要减去，

这和莫比乌斯函数的定义正好相反，所以可以预处理出莫比乌斯函数，这样统计以后最后的结果就是要求的答案。

所以就是求Σd * -mobius(d) * π(Ai/d)，但是现在复杂度是o(n^2)，所以要进行优化

考虑到π(Ai/d)中可能会有很多重复的数(例如，d = 5 时， Ai = 5, 6, 7, 8, 9时 Ai/d是相等的)，可以想办法将相同的数的个数统计出来，

然后使用快速幂计算多个数的乘积降低复杂度。所以从枚举最大公约数改为枚举 最大公约数的倍数 的个数, 即Ai/d,

下限和上限分别为,1和max(A)/d, 求Σd * -mobius(d) * π(t^f(t,d)), t = Ai/d, f(t,d)为t的个数, n + n/2 + n/3 + ... + n / n = n * (1 + 1/2 + 1/3 + ... + 1/n) ≈ n*log(n),再乘上快速幂的复杂度，总的复杂度为n * (log(n)) ^ 2.

计算f(t,d)时，可以将数组每个数Ai的个数统计出来，然后求前缀和，记为cnt[Ai]，对于每个d和t，f(d,t) = cnt[d * t + d - 1] - cnt[d * t - 1].

代码：

# include <iostream># include <algorithm># include <cstdio># include <cstring>using namespace std;typedef long long ll;const int maxn = 1e5 + 5;const int mod = 1e9 + 7;int vis[maxn];int pri[maxn];int mu[maxn];int a[maxn];int cnt[maxn * 2];int n;int len;void init() {    memset(vis, 0, sizeof vis);    mu[1] = 1; len = 0;    for (int i = 2; i < maxn; ++i) {        if (!vis[i]) {            pri[len++] = i;            mu[i] = -1;        }        for (int j = 0; j < len && i * pri[j] < maxn; ++j) {            vis[i * pri[j]] = 1;            if (i % pri[j]) mu[i * pri[j]] = -mu[i];            else {                mu[i * pri[j]] = 0;                break;            }        }    }}int fast_pow(int x, int n) {    int r = 1;    while (n) {        if (n & 1) r = (ll)r * x % mod;        x = (ll)x * x % mod;        n >>= 1;    }    return r;}int main(void){    init();    int T, Case = 0; scanf("%d", &T);    while (T-- && scanf("%d", &n)) {        memset(cnt, 0, sizeof cnt);        int Min = 1e5 + 1, Max = -1;        for (int i = 0; i < n; ++i) {            scanf("%d", a + i); ++cnt[a[i]];            Min = min(Min, a[i]); Max = max(Max, a[i]);        }        for (int i = Min; i <= Max * 2; ++i) cnt[i] += cnt[i - 1];        int ans = 0;        for (int d = 2; d <= Min; ++d) {            int tmp = 1;            for (int i = 1; i <= Max / d; ++i) {                tmp = (ll)tmp * fast_pow(i, cnt[i * d + d - 1] - cnt[i * d - 1]) % mod;            }            ans = ((ll)ans - mu[d] * tmp + mod) % mod;        }        printf("Case #%d: %d\n", ++Case, ans);    }    return 0;}