HDU 6053 容斥dp 或 莫比乌斯反演

题意：

题目链接：http://acm.hdu.edu.cn/showproblem.php?pid=6053
给出一段长度为n的序列a，要求设计同样长度为n的序列b，要求满足：
1. 1 <= b[i] <= a[i]
2. 对于任意区间[l,r]，都要满足gcd(a[l], a[l+1]…a[r])>=2
问有多少种可能的b数组

---

思路：

分析题意其实就是要求出有多少种方案可以使b数组全体gcd不为1。
有两种方式，可以利用莫比乌斯反演，也可以进行dp。

莫比乌斯反演：

设F[d]为n个数gcd为d的倍数的方案数，f[d]为n个数gcd==d的方案数。
只要求出F[d]就能计算f[d]，关键是F[d]的计算，枚举gcd==d，若对每个d都计算所有a[i]/d的乘积，显然会T。所以这里用了个技巧，只要知道a[i]/d相同的有多少个就行，这可以枚举d的倍数，然后里用前缀和求快速幂即可。

dp：

思路莫比乌斯反演类似，也是枚举gcd，这里可以设dp[i]为gcd==i的方案数，利用求F(d)的方案先求dp[d]，再利用容斥思想，将多余的部分都减掉即可，注意要逆序减掉。

---

代码：

莫比乌斯反演：

#include <bits/stdc++.h>
using namespace std;
typedef long long LL;
const LL MOD = 1e9 +7;
const LL INF = 0x3f3f3f3f;
const int MAXN = 1e5 + 10;

int u[MAXN], prime[MAXN];
bool vis[MAXN];

void Mobius(int x) {
    int cnt = 0;
    memset(vis, false, sizeof(vis));
    u[1] = 1;
    for (int i = 2; i <= x; i++) {
        if (!vis[i]) {
            prime[++cnt] = i;
            u[i] = -1;
        }
        for (int j = 1; j <= cnt && prime[j] * i <= x; j++) {
            vis[prime[j] * i] = true;
            if (i % prime[j]) u[prime[j] * i] = -u[i];
            else {
                u[prime[j] * i] = 0;
                break;
            }
        }
    }
}

LL pow_mod(LL a, LL n) {
    LL res = 1;
    while (n) {
        if (n & 1) res = res * a % MOD;
        a = a * a % MOD;
        n >>= 1;
    }
    return res;
}

int sum[MAXN];

int main() {
    //freopen("in.txt", "r", stdin);
    Mobius(100000);
    int T, cs = 0;
    scanf("%d", &T);
    while (T--) {
        int n, x, Min = INF, Max = 0;
        memset(sum, 0, sizeof(sum));
        scanf("%d", &n);
        LL ss = 1;
        for (int i = 1; i <= n; i++) {
            scanf("%d", &x);
            sum[x]++;
            Min = min(Min, x);
            Max = max(Max, x);
            ss = ss *x % MOD;
        }
        sum[0] = 0;
        for (int i = 1; i <= Max; i++) {
            sum[i] += sum[i - 1];
        }
        LL ans = 0;
        for (int i = 1; i <= Min; i++) {
            LL tmp = 1;
            for (int j = 1; j <= Max / i; j++) {
                int l = j * i, r = min((j + 1) * i - 1, Max);
                int cnt = sum[r] - sum[l - 1];
                tmp = tmp * pow_mod(j, cnt) % MOD;
            }

            ans = ((ans + u[i] * tmp % MOD) % MOD + MOD) % MOD;
        }
        ans = ((ss - ans) % MOD + MOD) % MOD;
        printf("Case #%d: %I64d\n", ++cs, ans);
    }
    return 0;
}

容斥+dp：

#include <bits/stdc++.h>
using namespace std;
typedef long long LL;
const LL MOD = 1e9 +7;
const LL INF = 0x3f3f3f3f;
const int MAXN = 1e5 + 10;

int dp[MAXN];

LL pow_mod(LL a, LL n) {
    LL res = 1;
    while (n) {
        if (n & 1) res = res * a % MOD;
        a = a * a % MOD;
        n >>= 1;
    }
    return res;
}

int sum[MAXN];

int main() {
    //freopen("in.txt", "r", stdin);
    int T, cs = 0;
    scanf("%d", &T);
    while (T--) {
        int n, x, Min = INF, Max = 0;
        memset(sum, 0, sizeof(sum));
        scanf("%d", &n);
        LL ss = 1;
        for (int i = 1; i <= n; i++) {
            scanf("%d", &x);
            sum[x]++;
            Min = min(Min, x);
            Max = max(Max, x);
            ss = ss *x % MOD;
        }
        sum[0] = 0;
        for (int i = 1; i <= Max; i++) {
            sum[i] += sum[i - 1];
        }
        for (int i = 1; i <= Min; i++) {
            dp[i] = 1;
            for (int j = 1; j <= Max / i; j++) {
                int l = j * i, r = min((j + 1) * i - 1, Max);
                int cnt = sum[r] - sum[l - 1];
                dp[i] = dp[i] * pow_mod(j, cnt) % MOD;
            }
        }
        for (int i = Min; i >= 2; i--) {
            for (int j = i * 2; j <= Min; j += i) {
                dp[i] = ((dp[i]- dp[j]) % MOD + MOD) % MOD;
            }
        }
        LL ans = 0;
        for (int i = 2; i <= Min; i++) {
            ans = (ans + dp[i]) % MOD;
        }
        printf("Case #%d: %I64d\n", ++cs, ans);
    }
    return 0;
}