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;
}