码迷，mamicode.com
首页 > 其他好文 > 详细

Mobius 反演

时间：2016-08-21 19:47:28      阅读：139      评论：0      收藏：0      [点我收藏+]

标签：

上次看莫比乌斯繁衍反演是一个月前，讲道理没怎么看懂..

然后出去跪了二十天， 然后今天又开始看发现其实并不难理解
 
开个这个仅记录一下写过的题。

 

HAOI 2011 B
 
这应该是莫比乌斯反演的模板题，有很多题解，不多说。
CODE：
//HAOI 2011 B
//by Cydiater
//2016.7.25
#include <iostream>
#include <cstring>
#include <string>
#include <algorithm>
#include <cstdio>
#include <cstdlib>
#include <ctime>
#include <cmath>
#include <iomanip>
#include <string>
#include <queue>
#include <map>
using namespace std;
#define ll long long
#define up(i,j,n)       for(int i=j;i<=n;i++)
#define down(i,j,n)     for(int i=j;i>=n;i--)
#define FILE "b"
const int MAXN=5e5+5;
const int oo=0x3f3f3f3f;
inline int read(){
      char ch=getchar();int x=0,f=1;
      while(ch>‘9‘||ch<‘0‘){if(ch==‘-‘)f=-1;ch=getchar();}
      while(ch>=‘0‘&&ch<=‘9‘){x=x*10+ch-‘0‘;ch=getchar();}
      return x*f;
}
int prime[MAXN],mu[MAXN],cnt=0,T,a,b,c,d,sum[MAXN],k;
bool vis[MAXN];
namespace solution{
      void make_mu(){
            mu[1]=1;
            memset(vis,0,sizeof(vis));
            up(i,2,50000){
                  if(!vis[i]){prime[++cnt]=i;mu[i]=-1;}
                  for(int j=1;prime[j]*i<=50000&&j<=cnt;j++){
                        vis[prime[j]*i]=1;
                        if(i%prime[j]!=0)mu[i*prime[j]]=-mu[i];
                        else{
                              mu[i*prime[j]]=0;
                              break;
                        }
                  }
            }
            up(i,1,50000)sum[i]=sum[i-1]+mu[i];
      }
      int get(int n,int m){
            if(n>m)swap(n,m);
            int pos,ans=0;
            for(int i=1;i<=n;i=pos+1){
                  pos=min(n/(n/i),m/(m/i));
                  ans+=(sum[pos]-sum[i-1])*(n/i)*(m/i);
            }
            return ans;
      }
}
int main(){
      //freopen("input.in","r",stdin);
      freopen(FILE".in","r",stdin);
      freopen(FILE".out","w",stdout);
      T=read();
      using namespace solution;
      memset(sum,0,sizeof(sum));
      make_mu();
      while(T--){
            a=read();b=read();c=read();d=read();k=read();
            a--;c--;
            a/=k;b/=k;c/=k;d/=k;
            printf("%d\n",get(b,d)-get(a,d)-get(b,c)+get(a,c));
      }
      return 0;
}

BZOJ 1101

上一道题的弱化版..
//BZOJ 1101
//by Cydiater
//2016.7.25
#include <iostream>
#include <cstdio>
#include <cstring>
#include <string>
#include <algorithm>
#include <ctime>
#include <cmath>
#include <cstdlib>
#include <queue>
#include <map>
#include <iomanip>
using namespace std;
#define ll long long
#define up(i,j,n)       for(int i=j;i<=n;i++)
#define down(i,j,n)     for(int i=j;i>=n;i--)
const int MAXN=5e4+50;
const int oo=0x3f3f3f3f;
inline int read(){
      char ch=getchar();int x=0,f=1;
      while(ch>‘9‘||ch<‘0‘){if(ch==‘-‘)f=-1;ch=getchar();}
      while(ch>=‘0‘&&ch<=‘9‘){x=x*10+ch-‘0‘;ch=getchar();}
      return x*f;
}
int mu[MAXN],prime[MAXN],cnt=0,sum[MAXN],T;
bool vis[MAXN];
namespace solution{
      void make_mu(){
            memset(mu,0,sizeof(mu));
            memset(sum,0,sizeof(sum));
            mu[1]=1;
            up(i,2,50000){
                  if(!vis[i]){prime[++cnt]=i;mu[i]=-1;}
                  for(int j=1;j<=cnt&&i*prime[j]<=50000;j++){
                        vis[i*prime[j]]=1;
                        if(i%prime[j]!=0)mu[i*prime[j]]=-mu[i];
                        else{
                              mu[i*prime[j]]=0;
                              break;
                        }
                  }
            }
            up(i,1,50000)sum[i]=sum[i-1]+mu[i];
      }
      int get(int n,int m){
            int pos,ans=0;
            if(n>m)swap(n,m);
            for(int i=1;i<=n;i=pos+1){
                  pos=min(m/(m/i),n/(n/i));
                  ans+=(sum[pos]-sum[i-1])*(m/i)*(n/i);
            }
            return ans;
      }
}
int main(){
      //freopen("input.in","r",stdin);
      using namespace solution;
      make_mu();
      T=read();
      while(T--){
            int a=read(),b=read(),k=read();
            a/=k;b/=k;
            printf("%d\n",get(a,b));
      }
      return 0;
}
Vijos 1889
 
题意就是让你求出第N个mu[i]!=0的数。
然后我们知道小于 i的莫比乌斯函数值不为0的数有Σmu[i]*(x/i^2)个（我从黄学长博客上看到的QAQ）
然后我们就可以愉快的二分了
 
PS注意枚举i时要设为ll
//Vijos 1889
//by Cydiater
//2016.7.25
#include <iostream>
#include <iomanip>
#include <cstdio>
#include <cstring>
#include <string>
#include <ctime>
#include <cmath>
#include <queue>
#include <map>
#include <cstdlib>
#include <algorithm>
using namespace std;
#define ll long long
#define up(i,j,n)       for(int i=j;i<=n;i++)
#define down(i,j,n)     for(int i=j;i>=n;i--)
const int MAXN=200000+50;
const int oo=0x3f3f3f3f;
inline ll read(){
      char ch=getchar();ll x=0,f=1;
      while(ch>‘9‘||ch<‘0‘){if(ch==‘-‘)f=-1;ch=getchar();}
      while(ch>=‘0‘&&ch<=‘9‘){x=x*10+ch-‘0‘;ch=getchar();}
      return x*f;
}
ll N,prime[MAXN],cnt=0,mu[MAXN],leftt,rightt,mid;
bool vis[MAXN];
namespace solution{
      void make_mu(){
            memset(vis,0,sizeof(vis));
            memset(mu,0,sizeof(mu));
            up(i,2,200000){
                  if(!vis[i]){prime[++cnt]=i;mu[i]=1;}
                  for(int j=1;prime[j]*i<=200000&&j<=cnt;j++){
                        vis[prime[j]*i]=1;
                        if(i%prime[j]!=0)mu[i*prime[j]]=-mu[i];
                        else{
                              mu[i*prime[j]]=0;
                              break;
                        }
                  }
            }
      }
      void init(){
            N=read();
      }
      bool check(ll num){
            ll lim=1LL*sqrt(1.0*num),ans=0;
            for(ll i=2;i<=lim;i++)
            ans+=num/(i*i)*mu[i];
            return ans>=N;
      }
      void slove(){
            leftt=N;rightt=25505460948LL;
            while(leftt+1<rightt){
                  mid=(leftt+rightt)>>1;
                  if(check(mid))    rightt=mid;
                  else              leftt=mid;
            }
            if(check(leftt))cout<<leftt<<endl;
            else            cout<<rightt<<endl;
      }
}
int main(){
      //freopen("input.in","r",stdin);
      //freopen("out.out","w",stdout);
      using namespace solution;
      make_mu();
      init();
      slove();
      return 0;
}

BZOJ 2440

上面那道题做一下小修改就好了
//bzoj2440
//by Cydiater
//2016.7.25
#include <iostream>
#include <iomanip>
#include <cstdio>
#include <cstring>
#include <string>
#include <ctime>
#include <cmath>
#include <queue>
#include <map>
#include <cstdlib>
#include <algorithm>
using namespace std;
#define ll long long
#define up(i,j,n)       for(int i=j;i<=n;i++)
#define down(i,j,n)     for(int i=j;i>=n;i--)
const int MAXN=200000+50;
const int oo=0x3f3f3f3f;
inline ll read(){
      char ch=getchar();ll x=0,f=1;
      while(ch>‘9‘||ch<‘0‘){if(ch==‘-‘)f=-1;ch=getchar();}
      while(ch>=‘0‘&&ch<=‘9‘){x=x*10+ch-‘0‘;ch=getchar();}
      return x*f;
}
ll N,prime[MAXN],cnt=0,mu[MAXN],leftt,rightt,mid,T;
bool vis[MAXN];
namespace solution{
      void make_mu(){
            memset(vis,0,sizeof(vis));
            memset(mu,0,sizeof(mu));
            up(i,2,200000){
                  if(!vis[i]){prime[++cnt]=i;mu[i]=1;}
                  for(int j=1;prime[j]*i<=200000&&j<=cnt;j++){
                        vis[prime[j]*i]=1;
                        if(i%prime[j]!=0)mu[i*prime[j]]=-mu[i];
                        else{
                              mu[i*prime[j]]=0;
                              break;
                        }
                  }
            }
      }
      void init(){
            N=read();
      }
      bool check(ll num){
            ll lim=1LL*sqrt(1.0*num),ans=0;
            for(ll i=2;i<=lim;i++)
            ans+=num/(i*i)*mu[i];
            return num-ans>=N;
      }
      void slove(){
            leftt=N;rightt=25505460948LL;
            while(leftt+1<rightt){
                  mid=(leftt+rightt)>>1;
                  if(check(mid))    rightt=mid;
                  else              leftt=mid;
            }
            if(check(leftt))cout<<leftt<<endl;
            else            cout<<rightt<<endl;
      }
}
int main(){
      //freopen("input.in","r",stdin);
      //freopen("out.out","w",stdout);
      using namespace solution;
      make_mu();
      T=read();
      while(T--){
            init();
            slove();
      }
      return 0;
}

YY的GCD

要克服公式恐惧症啊= =
 
和B那道题很像，但是与那道题不同的是这道题要求求出所有的素数。
设f(p)为在x在1-n中，y在1-m中满足gcd(x,y)==p的个数
所以我们就可以推出这样一个式子：
技术分享
 
然后我们就有了暴力求出这个题答案的方法了。
但是显然会超时，考虑优化
技术分享
 
 
这样的话如果我们能预处理出技术分享，就能很快的求出答案了。显然暴力的方法能预处理出来。
//YY de GCD
//by Cydiater
//2016.7.26
#include <iostream>
#include <cstring>
#include <string>
#include <algorithm>
#include <cmath>
#include <queue>
#include <map>
#include <cstdlib>
#include <cstdio>
#include <iomanip>
#include <ctime>
using namespace std;
#define ll long long
#define up(i,j,n)       for(int i=j;i<=n;i++)
#define down(i,j,n)     for(int i=j;i>=n;i--)
#define FILE "YYnoGCD"
const int oo=0x3f3f3f3f;
const int MAXN=1e7+5;
inline ll read(){
      char ch=getchar();ll x=0,f=1;
      while(ch>‘9‘||ch<‘0‘){if(ch==‘-‘)f=-1;ch=getchar();}
      while(ch>=‘0‘&&ch<=‘9‘){x=x*10+ch-‘0‘;ch=getchar();}
      return x*f;
}
int prime[MAXN],cnt=0,f[MAXN],T,mu[MAXN];
bool vis[MAXN];
namespace solution{
      void make_mu(){
            memset(vis,0,sizeof(vis));
            memset(f,0,sizeof(f));
            mu[1]=1;
            up(i,2,10000000){
                  if(!vis[i]){prime[++cnt]=i;mu[i]=-1;}
                  up(j,1,cnt){
                        if(prime[j]*i>10000000)break;
                        vis[prime[j]*i]=1;
                        if(i%prime[j]!=0)mu[i*prime[j]]=-mu[i];
                        else{
                              mu[i*prime[j]]=0;
                              break;
                        }
                  }
            }
            up(i,1,cnt)up(j,1,10000000){
                  if(prime[i]*j>10000000)break;
                  f[prime[i]*j]+=mu[j];
            }
            up(i,1,10000000)f[i]+=f[i-1];
      }
      void slove(ll a,ll b){
            ll ans=0,pos;
            if(a>b)swap(a,b);
            up(i,1,a){
                  pos=min(a/(a/i),b/(b/i));
                  ans+=(f[pos]-f[i-1])*(a/i)*(b/i);
                  i=pos;
            }
            cout<<ans<<endl;
      }
}
int main(){
      //freopen("input.in","r",stdin);
      freopen(FILE".in","r",stdin);
      freopen(FILE".out","w",stdout);
      using namespace solution;
      make_mu();
      T=read();
      while(T--){
            ll a=read(),b=read();
            slove(a,b);
      }
      return 0;
}

 

BZOJ 4407于神之怒

 
这些题为什么越来越难了QAQ
给下N,M,K.求
技术分享
下面给出公式的推导：
技术分享
 
技术分享
  
//BZOJ 4407
//by Cydiater
//2016.7.27
#include <iostream>
#include <cstdio>
#include <cstring>
#include <string>
#include <algorithm>
#include <cstdlib>
#include <queue>
#include <map>
#include <iomanip>
#include <cmath>
#include <ctime>
using namespace std;
#define ll long long
#define up(i,j,n)       for(int i=j;i<=n;i++)
#define down(i,j,n)     for(int i=j;i>=n;i--)
const int MAXN=5e6+5;
const int LIM=5e6;
const int oo=0x3f3f3f3f;
const int mod=1e9+7;
inline int read(){
      char ch=getchar();int x=0,f=1;
      while(ch>‘9‘||ch<‘0‘){if(ch==‘-‘)f=-1;ch=getchar();}
      while(ch>=‘0‘&&ch<=‘9‘){x=x*10+ch-‘0‘;ch=getchar();}
      return x*f;
}
ll f[MAXN],prime[MAXN],cnt=0,N,K,g[MAXN],T;
bool vis[MAXN];
namespace solution{
      inline ll quick_pow(ll a,ll b){
            ll tmp=1;
            while(b){
                  if(b&1)tmp=(tmp*a)%mod;
                  a=(a*a)%mod;
                  b>>=1;
            }
            return tmp;
      }
      void pret(){
            memset(vis,0,sizeof(vis));
            f[1]=1;
            up(i,2,LIM){
                  if(!vis[i]){prime[++cnt]=i;g[cnt]=quick_pow(i,K);f[i]=g[cnt]-1;}
                  up(j,1,cnt){
                        if(i*prime[j]>LIM)break;
                        vis[i*prime[j]]=1;
                        if(i%prime[j]!=0){
                              f[i*prime[j]]=(f[i]*f[prime[j]])%mod;
                        }else{
                              f[i*prime[j]]=(f[i]*g[j])%mod;
                              break;
                        }
                  }
            }
            up(i,1,LIM)f[i]=(f[i-1]+f[i])%mod;
      }
      void slove(ll a,ll b){
            ll ans=0,pos;
            if(a>b)swap(a,b);
            up(i,1,a){
                  pos=min(a/(a/i),b/(b/i));
                  ans+=(((f[pos]+mod-f[i-1])*(a/i))%mod)*(b/i);
                  ans%=mod;
                  i=pos;
            }
            printf("%lld\n",ans);
      }
}
int main(){
      freopen("input.in","r",stdin);
      using namespace solution;
      T=read();K=read();
      pret();
      while(T--)slove(read(),read());
      return 0;
}


BZOJ 2154

数论好坑啊
做这道题是为了做下一道题，这道题是下一道题的弱化版..但是我调试了几乎一个上午。
给定N,M求N,M内的lcm累加和。
下面是公式时间：
技术分享
 
 

 技术分享
 然后双重分块就好了

这个调试的真恶心...
  
//BZOJ 2154
//by Cydiater
//2016.7.27
#include <iostream>
#include <cstdio>
#include <string>
#include <cstring>
#include <iomanip>
#include <queue>
#include <map>
#include <ctime>
#include <cmath>
#include <algorithm>
#include <cstdlib>
using namespace std;
#define ll long long
#define up(i,j,n)       for(ll i=j;i<=n;i++)
#define down(i,j,n)     for(ll i=j;i>=n;i--)
const ll MAXN=1e7+5;
const ll LIM=1e7;
const ll mod=20101009;
const ll oo=0x3f3f3f3f;
inline ll read(){
      char ch=getchar();ll x=0,f=1;
      while(ch>‘9‘||ch<‘0‘){if(ch==‘-‘)f=-1;ch=getchar();}
      while(ch>=‘0‘&&ch<=‘9‘){x=x*10+ch-‘0‘;ch=getchar();}
      return x*f;
}
ll prime[MAXN],cnt=0,mu[MAXN],N,M;
bool vis[MAXN];
namespace solution{
      void pret(){
            memset(vis,0,sizeof(vis));
            memset(mu,0,sizeof(mu));
            mu[1]=1;
            up(i,2,N){
                  if(!vis[i]){prime[++cnt]=i;mu[i]=-1;}
                  up(j,1,cnt){
                        if(prime[j]*i>N)break;
                        vis[prime[j]*i]=1;
                        if(i%prime[j]!=0){
                              mu[i*prime[j]]=-mu[i];
                        }else{
                              mu[i*prime[j]]=0;
                              break;
                        }
                  }
            }
            up(i,1,N)mu[i]=(mu[i-1]+(mu[i]*i*i)%mod)%mod;
      }
      ll sum(ll a,ll b){
            return ((a+1)*a/2%mod)*(b*(b+1)/2%mod)%mod;
      }
      ll F(ll a,ll b){
            ll pos,ans=0;
            if(a>b)swap(a,b);
            up(i,1,a){
                  pos=min(a/(a/i),b/(b/i));
                  ans=(ans+(mu[pos]-mu[i-1])*sum(a/i,b/i))%mod;
                  i=pos;
            }
            return ans;
      }
      void slove(ll a,ll b){
            if(a>b)swap(a,b);
            ll pos,ans=0;
            up(i,1,a){
                  pos=min(a/(a/i),b/(b/i));
                  ans=(ans+(i+pos)*(pos-i+1)/2%mod*F(a/i,b/i)%mod)%mod;
                  i=pos;
            }
            printf("%lld\n",(ans+mod)%mod);
      }
}
int main(){
      freopen("input.in","r",stdin);
      using namespace solution;
      N=read();M=read();
      if(N>M)swap(N,M);
      pret();
      slove(N,M);
      //cout<<"Time has passed:"<<1.0*clock()/1000<<"s!"<<endl;
      return 0;
}

BZOJ 2693:

这些数论题套路很多啊..
接着上一个公式：
技术分享
//BZOJ 2693
//by Cydiater
//2016.7.29
#include <iostream>
#include <cstdio>
#include <cstring>
#include <string>
#include <algorithm>
#include <cmath>
#include <ctime>
#include <queue>
#include <map>
#include <cstdlib>
#include <iomanip>
using namespace std;
#define ll long long
#define up(i,j,n)       for(ll i=j;i<=n;i++)
#define down(i,j,n)     for(ll i=j;i>=n;i--)
const int MAXN=1e7+10;
const int LIM=1e7;
const int mod=1e8+9;
const int oo=0x3f3f3f3f;
inline ll read(){
      char ch=getchar();ll x=0,f=1;
      while(ch>‘9‘||ch<‘0‘){if(ch==‘-‘)f=-1;ch=getchar();}
      while(ch>=‘0‘&&ch<=‘9‘){x=x*10+ch-‘0‘;ch=getchar();}
      return x*f;
}
ll f[MAXN],prime[MAXN],cnt=0,T;
bool vis[MAXN];
namespace solution{
      void pret(){
            memset(vis,0,sizeof(vis));
            f[1]=1;
            up(i,2,LIM){
                  if(!vis[i]){prime[++cnt]=i;f[i]=(i-i*i%mod)%mod;}
                  up(j,1,cnt){
                        if(prime[j]*i>LIM)break;
                        vis[prime[j]*i]=1;
                        if(i%prime[j]==0){
                              f[i*prime[j]]=(f[i]*prime[j])%mod;
                              break;
                        }
                        f[i*prime[j]]=(f[prime[j]]*f[i])%mod;
                  }
            }
            up(i,1,LIM)f[i]=(f[i]+f[i-1])%mod;
      }
      ll sum(ll a,ll b){
            ll tmp1=(a*(a+1)/2)%mod;
            ll tmp2=(b*(b+1)/2)%mod;
            return (tmp1*tmp2)%mod;
      }
      void slove(ll a,ll b){
            if(a>b)swap(a,b);
            ll pos,ans=0;
            up(i,1,a){
                  pos=min(a/(a/i),b/(b/i));
                  ans=(ans+(f[pos]-f[i-1])*sum(a/i,b/i)%mod)%mod;
                  i=pos;
            }
            printf("%lld\n",(ans+mod)%mod);
      }
}
int main(){
      //freopen("input.in","r",stdin);
      using namespace solution;
      T=read();
      pret();
      while(T--)slove(read(),read());
      return 0;
}

 

Mobius 反演

标签：

原文地址：http://www.cnblogs.com/Cydiater/p/5793255.html
踩
(0)
赞
(0)
   
举报
评论 一句话评论（0）
登录后才能评论！
分享档案
更多>
2021年07月29日 (22)
2021年07月28日 (40)
2021年07月27日 (32)
2021年07月26日 (79)
2021年07月23日 (29)
2021年07月22日 (30)
2021年07月21日 (42)
2021年07月20日 (16)
2021年07月19日 (90)
2021年07月16日 (35)
周排行
mamicode.com排行更多图片
更多
友情链接
兰亭集智   国之画   百度统计   站长统计   阿里云   chrome插件   新版天听网
关于我们 - 联系我们 - 留言反馈
© 2014 mamicode.com 版权所有  联系我们:gaon5@hotmail.com
迷上了代码！