hdu 6053 TrickGCD

Problem

acm.hdu.edu.cn/showproblem.php?pid=6053

莫比乌斯反演介绍

莫比乌斯反演（by PoPoQQQ）
莫比乌斯反演入门
莫比乌斯反演

Reference

2017 Multi-University Training Contest 2 solutions BY 电子科技大学（1009题）
HDU 6053 TrickGCD （莫比乌斯函数）
解题报告：HDU_6053 TrickGCD 莫比乌斯反演

Meaning

给定一个序列 {An}，要构造一个序列 {Bn}，满足：

• 1≤Bi≤Ai
• gcd（B1，…，Bn）≥2

问构造的方案数。

Analysis

暴力的做法，就是枚举 gcd，每一个 gcd 产生的答案是 ∏ni=1Aigcd，总的答案就是它们的求和，但要去重，可以用莫比乌斯反演。
…好吧其实现在并不太懂莫比乌斯反演…按官方题解，定义：
F(x)：gcd 是 x 倍数的方案数；
我猜相应的就有：
f(x)：gcd 是 x 的方案数；
于是：
F(x)=∑x|df(d)⇒f(x)=∑x|dμ(dx)F(d)
而 F(x)=∏ni=1Aix，就可以做了。
按 dalao 的说法，因为枚举 gcd 时是从 2 开始的，f(1) 没有用到，所以在乘 μ(i)的时候，要多乘一个 -1，于是就有了第二个参考博客的反演式。
在处理那个连乘式的时候，要用一点优化处理。
好吧还是很迷…

Code

#include <cstdio>#include <cstring>#include <algorithm>using namespace std;const int N = 100000, A = N, MOD = 1000000007;int a[N], sum[N<<1|1];int prime[A+1];int mu[A+1];void Mobius(){    mu[1] = 1;    memset(prime, -1, sizeof prime);    for(int i = 2, top = 0; i <= A; ++i)    {        if(prime[i])        {            prime[top++] = i;            mu[i] = -1;        }        for(int j = 0; j < top && i * prime[j] <= A; ++j)        {            prime[i * prime[j]] = 0;            if(i % prime[j])                mu[i * prime[j]] = -mu[i];            else            {                mu[i * prime[j]] = 0;                break;            }        }    }}int fast_pow(long long a, int b){    long long x = 1;    for( ; b; b >>= 1, a = a * a % MOD)        if(b & 1)            x = x * a % MOD;    return x % MOD;}int main(){    Mobius();    int T;    scanf("%d", &T);    for(int kase = 1; kase <= T; ++kase)    {        int n, small = A, big = 0;        scanf("%d", &n);        memset(sum, 0, sizeof sum);        for(int i = 0; i < n; ++i)        {            scanf("%d", a+i);            ++sum[a[i]];            small = min(small, a[i]);            big = max(big, a[i]);        }        for(int i = 1; i <= big << 1; ++i)            sum[i] += sum[i-1];        long long ans = 0, tmp;        for(int d = 2; d <= small; ++d)            if(tmp = -mu[d]) // mu[d] != 0            {                for(int i = d; i <= big; i += d)                    tmp = tmp * fast_pow(i/d, sum[i+d-1]-sum[i-1]) % MOD;                ans = (ans + tmp + MOD) % MOD;            }        printf("Case #%d: %I64d\n", kase, ans);    }    return 0;}