PoPoQQQ大神的PPT:
https://wenku.baidu.com/view/fbec9c63ba1aa8114431d9ac.html?from=search
其实莫比乌斯反演就这个东西:
如果有:
那么有:
另一种形式是如果
那么
其中
(1)
μ(1)=1
(2)若d为若干个互异素数之积,即
d=p1p2p3...pk
,那么有
μ(d)=(−1)k
(3)其他情况下
μ(d)=0
μ
是积性函数,即对于任意正整数x,y,如果有gcd(x,y)=1,那么有
μ(x)μ(y)=μ(xy)
所以我们可以用线性筛求出
μ
:
mu[1]=1;
for(int i=2;i<=x;i++){
if(!ok[i]){
prime[++prime[0]]=i;
mu[i]=-1;
}
for(int j=1;j<=prime[0]&&i*prime[j]<=x;j++){
ok[i*prime[j]]=true;
if(i%prime[j])mu[i*prime[j]]=-mu[i];
else{
mu[i*prime[j]]=0;
break;
}
}
}
所以我们有了它能做什么?
如果有一个不是很好求的函数f,那么我们可能可以通过反演求出它的值.
首先一道最简单的例题:Bzoj1101 Zap
这是个权限题,看不了没关系,题目大意就是对于给定的整数a,b和d,有多少正整数对x,y,满足
x≤a,y≤b
,并且gcd(x,y)=d。
即求:
我们可以把d除掉,令 n=⌊ad⌋,m=⌊bd⌋ ,那么答案就是:
我们不妨令
然后我们会发现,
莫比乌斯反演得
而F(i)是很好求的,因为有
所以我们要的答案f(1)就是等于:
现在我们可以O(n)出解,但是由于多组数据我们还要以优化.
考虑到 ⌊ni⌋⌊mi⌋ 最多只有 2(n√+m‾‾√) 个取值,所以我们按照这个的值分块计算,所以预处理一下 μ 的前缀和即可.
那么这个代码怎么写?
LL nxt,tot=0;
for(int i=1;i<=min(n,m);i=nxt+1){
nxt=min(n/(n/i),m/(m/i));//计算当前i对应的值域的上界
tot+=(n/i)*(m/i)*(sum[nxt]-sum[i-1]);
//sum是莫比乌斯函数的前缀和
}
return tot;
所以现在复杂度是O(
n√
),多组数据可以过.
不过问题来了,我们是怎么想到F这个函数对进行反演的?如果在题目简单的情况下,可以比较容易看出F,可如果题目难,我们怎么知道用哪个函数来反演?
我们可以先试着用容斥的方法做这道题:
范围内gcd为1的数对个数=(范围内gcd是1的倍数的数对个数)-(范围内gcd是2的倍数的数对个数)-(范围内gcd是3的倍数的数对个数)-(范围内gcd是5的倍数的数对个数)+(范围内gcd是6的倍数的数对个数)……
因为(范围内gcd是4的倍数的数对个数)在第2项处已经被减掉了,所以不用再减一次.
而(范围内gcd是6的倍数的数对个数)在第2项和第3项处被减掉了2次,所以要加回来一次.
我们发现了什么?是不是这里有一个莫比乌斯函数?
是不是很神奇?莫比乌斯函数就是容斥中每一项的系数!
其实根据莫比乌斯函数的定义和容斥的规则我们就可以发现这个性质.
所以这题我们可以根据容斥直接列式计算,上面的容斥就是这个式子:
同样,很多莫比乌斯反演的题目都可以根据容斥直接列式.
Bzoj2301 Problem B
这题比上一题多了一个下界,我们转化一下答案即可用同样的方法解决.
代码:
#include<cstdio>
#include<algorithm>
typedef long long LL;
const int maxn=50010;
int T,a,b,c,d,k,mu[maxn],prime[maxn],sum[maxn];
bool not_prime[maxn];
void Make(){
not_prime[1]=mu[1]=1;
for(int i=1;i<=50000;i++){
if(!not_prime[i]){
prime[++prime[0]]=i;
mu[i]=-1;
}
for(int j=1;j<=prime[0]&&i*prime[j]<=50000;j++){
not_prime[i*prime[j]]=true;
if(i%prime[j]==0){
mu[i*prime[j]]=0;
break;
}
mu[i*prime[j]]=-mu[i];
}
}
for(int i=1;i<=50000;i++)sum[i]=sum[i-1]+mu[i];
}
LL Work(LL n,LL m){
n/=k;m/=k;
if((!n)||(!m))return 0;
if(n>m)std::swap(n,m);
LL tot=0;
for(int i=1,nxt;i<=n;i=nxt+1){
nxt=std::min(n/(n/i),m/(m/i));
tot+=(n/i)*(m/i)*(sum[nxt]-sum[i-1]);
}
return tot;
}
int main(){
Make();
scanf("%d",&T);
while(T--){
scanf("%d%d%d%d%d",&a,&b,&c,&d,&k);
printf("%lld\n",Work(b,d)-Work(a-1,d)-Work(b,c-1)+Work(a-1,c-1));
}
return 0;
}
Bzoj2440 完全平方数
考虑二分答案.对于当前答案n,我们考虑如何求出它之前(包括它)有多少个数是合法的.用容斥列式,枚举完全平方数:
其中 ⌊ni2⌋ 是n以内是 i2 整数倍的数的个数.
由于 i2 大于n时 ⌊ni2⌋ 为0,所以我们只需要枚举i到 n√ 即可.
复杂度O( n√log2n )
代码:
#include<cstdio>
#include<algorithm>
using namespace std;
typedef long long LL;
const int maxn=100010,maxl=100000;
int T,n,mu[maxn],prime[maxn];
bool ok[maxn];
void Make(){
mu[1]=1;
for(int i=2;i<=maxl;i++){
if(!ok[i]){
prime[++prime[0]]=i;
mu[i]=-1;
}
for(int j=1;j<=prime[0]&&i*prime[j]<=maxl;j++){
ok[i*prime[j]]=true;
if(i%prime[j])mu[i*prime[j]]=-mu[i];
else{
mu[i*prime[j]]=0;
break;
}
}
}
}
LL check(LL ans){
LL tot=0;
for(LL i=1;i*i<=ans;i++){
tot+=ans/(i*i)*mu[i];
}
return tot;
}
LL Binary(LL l,LL r){
LL mid,ans;
while(l<=r){
mid=(l+r)>>1;
if(check(mid)<n)l=mid+1;
else ans=mid,r=mid-1;
}
return ans;
}
int main(){
Make();
scanf("%d",&T);
while(T--){
scanf("%d",&n);
printf("%lld\n",Binary(n,n<<1));
}
return 0;
}
Bzoj3529 数表
大意:
F(i)为i的约数和,求:
要求计算进入结果的 F(gcd(i,j))≤a
我们先不考虑a的限制.
令g(i)为
1≤x≤n,1≤y≤m
且gcd(x,y)=i的数对(x,y)的个数.
前面我们知道了
g(i)=∑i|dμ(di)⌊nd⌋⌊md⌋
,
则:
我们对于每个d预处理出 ∑i|dF(i)μ(di) 及其前缀和,这个可以O( nlog2n )做到,然后单组询问O( n√ ),可以过.
再考虑加了a的限制怎么做.
我们知道对当前答案有贡献的都是满足
所以不妨离线所有询问,按照a的增序回答询问,用树状数组维护 ∑i|dF(i)μ(di) 的前缀和,询问的a变大就新插入一些对答案有贡献的值.总复杂度O( nlog2n+qn√log2n ).
代码:
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
inline int read(){
int x=0,f=1;char ch=getchar();
for(;ch<'0'||ch>'9';ch=getchar())if(ch=='-')f=-1;
for(;ch>='0'&&ch<='9';ch=getchar())x=x*10+ch-'0';
return x*f;
}
typedef long long LL;
const int maxn=100010,maxl=100000;
const LL mod=1ll<<31;
struct array{
#define lowbit(x) (x&-x)
int n;
LL a[maxn];
array(){memset(a,0,sizeof a);n=maxl;}
void add(int x,LL k){
while(x<=n){
a[x]+=k;
x+=lowbit(x);
}
}
LL sum(int x){
LL tot=0;
while(x){
tot+=a[x];
x-=lowbit(x);
}
return tot;
}
}a;
struct Q{
LL n,m,a,id;
bool operator<(Q b)const{
return a<b.a;
}
}q[maxn];
int T,mu[maxn],prime[maxn];
LL ans[maxn];
pair<LL,int>f[maxn];
bool ok[maxn];
void Make(){
mu[1]=1;
for(int i=2;i<=maxl;i++){
if(!ok[i]){
prime[++prime[0]]=i;
mu[i]=-1;
}
for(int j=1;j<=prime[0]&&i*prime[j]<=maxl;j++){
ok[i*prime[j]]=true;
if(i%prime[j])mu[i*prime[j]]=-mu[i];
else{
mu[i*prime[j]]=0;
break;
}
}
}
for(int i=1;i<=maxl;i++){
for(int j=i;j<=maxl;j+=i){
f[j].first+=i;
}
}
for(int i=1;i<=maxl;i++)f[i].second=i;
}
LL Query(LL n,LL m){
if(n>m)swap(n,m);
LL tot=0;
for(int i=1,nxt;i<=n;i=nxt+1){
nxt=min(n/(n/i),m/(m/i));
tot+=(n/i)*(m/i)*(a.sum(nxt)-a.sum(i-1))%mod;
tot%=mod;tot=(tot+mod)%mod;
}
return (tot+mod)%mod;
}
int main(){
Make();T=read();
for(int i=1;i<=T;i++){
q[i].n=read();q[i].m=read();
q[i].a=read();q[i].id=i;
}
sort(f+1,f+maxl+1);sort(q+1,q+T+1);
for(int t=1,i=1;t<=T;t++){
for(;i<=maxl&&f[i].first<=q[t].a;i++){
for(int j=f[i].second;j<=maxl;j+=f[i].second){
a.add(j,mu[j/f[i].second]*f[i].first);
}
}
ans[q[t].id]=Query(q[t].n,q[t].m);
}
for(int i=1;i<=T;i++)printf("%lld\n",ans[i]);
return 0;
}
Bzoj2154 Crash的数字表格
大意:给定n,m(n,m
变形:
令
枚举d=gcd(i,j):
根据容斥可以列出F的表达式:
预处理出 μ(i)∗i2 的前缀和可以O( n√ )求出F.
所以可以O( n√n√ )=O(n)求出答案.
代码:
#include<cstdio>
#include<algorithm>
using namespace std;
typedef long long LL;
const int maxn=10000010,mod=20101009;
LL n,m,fac[maxn];
int mu[maxn],prime[maxn];
bool ok[maxn];
void Init(int x){
mu[1]=1;
for(int i=2;i<=x;i++){
if(!ok[i]){
prime[++prime[0]]=i;
mu[i]=-1;
}
for(int j=1;j<=prime[0]&&i*prime[j]<=x;j++){
ok[i*prime[j]]=true;
if(i%prime[j])mu[i*prime[j]]=-mu[i];
else{
mu[i*prime[j]]=0;
break;
}
}
}
for(LL i=1;i<=x;i++)fac[i]=((fac[i-1]+mu[i]*i*i%mod)%mod+mod)%mod;
}
LL F(LL x,LL y){
if(x>y)swap(x,y);
if(!x)return 0;
LL nxt,tot=0;
for(LL i=1;i<=x;i=nxt+1){
nxt=min(x/(x/i),y/(y/i));
tot=(tot+(x/i)*(x/i+1)/2%mod*((y/i)*(y/i+1)/2%mod)%mod*(((fac[nxt]-fac[i-1])%mod+mod)%mod)%mod)%mod;
}
return tot;
}
LL Calc(LL x,LL y){
LL nxt,tot=0;
for(LL i=1;i<=x;i=nxt+1){
nxt=min(x/(x/i),y/(y/i));
(tot+=(nxt+i)*(nxt-i+1)/2%mod*F(x/i,y/i)%mod)%=mod;
}
return tot;
}
int main(){
scanf("%lld%lld",&n,&m);
if(n>m)swap(n,m);
if(!n)return printf("0\n"),0;
Init(n);
printf("%lld\n",Calc(n,m));
return 0;
}
Bzoj2693 jzptab
上一题的多组数据版本,还要继续优化.
记
现在我们只要能O(n)之内求出 ∑i|Dμ(i)i 就可以O( n√ )完成一组询问.
令 f(i)=μ(i)i ,若x,y互质,那么
所以f为积性函数.又积性函数的约数和为积性函数,
所以 ∑i|Dμ(i)i 为积性函数,所以我们可以通过线性筛求出它的值.
此题到此解决,复杂度为O( n+qn√ ).
代码:
#include<cstdio>
#include<algorithm>
using namespace std;
typedef long long LL;
const int maxn=10000010,mod=100000009;
LL n,m,f[maxn],g[maxn],prime[maxn];
int T;
bool ok[maxn];
void Init(int x){
f[1]=1;
for(LL i=2;i<=x;i++){
if(!ok[i]){
prime[++prime[0]]=i;
f[i]=1-i+mod;
}
for(int j=1;j<=prime[0]&&i*prime[j]<=x;j++){
ok[i*prime[j]]=true;
if(i%prime[j])f[i*prime[j]]=(f[i]*(1-prime[j])%mod)+mod;
else{
f[i*prime[j]]=f[i];
break;
}
}
}
for(LL i=1;i<=x;i++)g[i]=(g[i-1]+i*f[i]%mod)%mod;
}
LL Sum(LL x,LL y){
return (x*(x+1)/2%mod)*(y*(y+1)/2%mod)%mod;
}
LL Calc(LL x,LL y){
LL nxt,tot=0;
for(LL i=1;i<=x;i=nxt+1){
nxt=min(x/(x/i),y/(y/i));
(tot+=Sum(x/i,y/i)*((((g[nxt]-g[i-1])%mod)+mod)%mod)%mod)%=mod;
}
return tot;
}
int main(){
Init(10000000);
scanf("%d",&T);
while(T--){
scanf("%lld%lld",&n,&m);
if(n>m)swap(n,m);
if(!n){
printf("0\n");
continue;
}
printf("%lld\n",Calc(n,m));
}
return 0;
}
Spoj LCM Sum
求:
这道题根上面两题看起来有点像,但是解法不同.我们先化简式子:
枚举d=gcd(i,n)
我们考虑怎么求出里面那个 ∑ 的值.
里面的那部分其实就是要求1~x内与x互质的数之和.
考虑反演,用容斥枚举gcd列式:
当x=1时整个式子等于1
当x>1时第一项为0,我们考虑第二项:
其实第二项中的 ∑d|xμ(d)d=ϕ(x)x
这个结论要记一下,证明如下:
所以两边同时除以x就得到 ∑d|xμ(d)d=ϕ(x)x
综上,当x>1时,1~x中与x互质的数的个数为
所以,
由于 nd=1 时 ∑ 里面的那一部分的值应该为1,但是我们直接算出来的值是 12 向下取整之后是0,所以我们可以人为的给 ∑ 里面加上一个1,即:
所以我们可以O(n)筛出欧拉函数,然后O( nlog2n )预处理出 ∑d|nndϕ(nd)2 ,就可以做到单组询问O(1).于是这道题就解决了.
代码:
#include<cstdio>
#include<algorithm>
using namespace std;
typedef long long LL;
const int maxn=1000010;
int T,phi[maxn],prime[maxn];
LL n,ans[maxn];
bool ok[maxn];
void Init(int x){
phi[1]=1;
for(int i=2;i<=x;i++){
if(!ok[i]){
prime[++prime[0]]=i;
phi[i]=i-1;
}
for(int j=1;j<=prime[0]&&i*prime[j]<=x;j++){
ok[i*prime[j]]=true;
if(i%prime[j])phi[i*prime[j]]=phi[i]*(prime[j]-1);
else{
phi[i*prime[j]]=phi[i]*prime[j];
break;
}
}
}
for(int i=1;i<=x;i++){
for(int j=i;j<=x;j+=i){
ans[j]+=(long long)phi[j/i]*(j/i)/2;
}
}
}
int main(){
Init(1000000);
scanf("%d",&T);
while(T--){
scanf("%lld",&n);
printf("%lld\n",n*ans[n]+n);
}
return 0;
}