莫比乌斯反演
基本内容:
线性求莫比乌斯函数模板:
#include<cstdio> #define N 10000000 using namespace std; int i,j,n,p,ans,R,cnt[N],prime[N],phi[N],u[N]; bool f[N]; int main(){ scanf("%d",&n);scanf("%d",&R);//R为逆元模 u[1]=1; for(i=2;i<=n;i++){ if(!f[i]){ prime[++p]=i;//prime存素数 cnt[i]=1;//cnt为因数个数 phi[i]=i-1;//phi为欧拉函数 u[i]=-1;//u为莫比乌斯函数 } for(j=1;j<=p&&i*prime[j]<=n;j++){ f[i*prime[j]]=1; cnt[i*prime[j]]=cnt[i]+1; if(i%prime[j]==0){ phi[i*prime[j]]=phi[i]*prime[j]; break; } phi[i*prime[j]]=phi[i]*(prime[j]-1); u[prime[j]*i]=-u[i]; } } }
BZOJ2301:莫比乌斯函数另一种形式+前缀和
i-(n/(n/i))可以取出和n/i商相同的值快速更新答案
#include<cstdio> #include<iostream> #define N 111111 using namespace std; long long i,j,T,a,b,c,d,k,last,p,u[N],sum[N],prime[N]; bool f[N]; long long calc(long long n,long long m){ n/=k;m/=k; long long re=0; for(long long i=1;i<=min(n,m);i=last+1){ last=min(n/(n/i),m/(m/i)); re+=(n/i)*(m/i)*(sum[last]-sum[i-1]); } return re; } int main(){ sum[1]=u[1]=1; for(i=2;i<N;i++){ if(!f[i]){ prime[++p]=i; u[i]=-1; } for(j=1;j<=p&&i*prime[j]<N;j++){ f[i*prime[j]]=1; if(i%prime[j])u[prime[j]*i]-=u[i]; else break; } sum[i]=sum[i-1]+u[i]; } scanf("%lld",&T); while(T--){ scanf("%lld%lld%lld%lld%lld",&a,&b,&c,&d,&k); printf("%lld\n",calc(b,d)+calc(a-1,c-1)-calc(a-1,d)-calc(b,c-1)); } }
#include<cstdio> #include<algorithm> #define N 50050 using namespace std; int n,m,i,j,pos,ans,t,Q,x,y,z,p[N],sum[N],u[N];bool f[N]; int main(){ for(n=50000,sum[1]=u[1]=1,i=2;i<=n;i++){ if(!f[i])p[++t]=i,u[i]=-1; for(j=1;j<=t&&i*p[j]<=n;j++){ f[i*p[j]]=1; if(i%p[j])u[i*p[j]]=-u[i];else break; } sum[i]=sum[i-1]+u[i]; } for(scanf("%d",&Q);Q--;printf("%d\n",ans),ans=0){ scanf("%d%d%d",&x,&y,&z);n=x/z;m=y/z;if(n>m)swap(n,m); for(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); } } }
BZOJ2820 枚举每个质数,那么答案和上题一样求,可是这样会T
使pd=T,那么ans=Σ(n/T)(m/T)Σu(T/p)(p|T),这样预处理使得时间复杂度变为O(n+qn^0.5)
#include<bits/stdc++.h> #define N 10001000 using namespace std; int i,j,t,T,n,m,la,p[N],sum[N],u[N];long long ans;bool f[N]; int main(){ for(u[1]=1,i=2;i<N;i++){ if(!f[i])p[++t]=i,u[i]=-1; for(j=1;j<=t&&i*p[j]<N;j++){ f[i*p[j]]=1; if(i%p[j]==0)break; u[i*p[j]]=-u[i]; } } for(j=1;j<=t;j++)for(i=p[j];i<N;i+=p[j])sum[i]+=u[i/p[j]]; for(i=1;i<N;i++)sum[i]+=sum[i-1]; for(scanf("%d",&T);T--;printf("%lld\n",ans),ans=0) for(scanf("%d%d",&n,&m),i=1;i<=min(n,m);i=la+1) ans+=(long long)(n/i)*(m/i)*(sum[la=min(n/(n/i),m/(m/i))]-sum[i-1]); }
BZOJ3529 用g(i)表示1≤x≤n,1≤y≤m,gcd(x,y)=i的个数,那么g(i)=Σu(d/i)*(n/d)*(m/d)
可得公式ans=ΣF(i)g(i)=(n/d)*(m/d)*ΣF(i)Σu(d/i),现在只需要求出ΣF(i)Σu(d/i)的前缀和即可
如果F(i)没有限制,那么可以O(nlgn)前缀和预处理,O(n^0.5)处理每个询问
现在有限制,可与对F(i)和所有a排序,对于每个询问a只把F(i)≤a的加入树状数组中,这样维护前缀和,总时间复杂度O(nlg^2n+qn^0.5lgn)
#include<bits/stdc++.h> #define N 100100 using namespace std; int i,j,k,t,T,u[N],p[N],c[N],ans[N];bool f[N]; struct P{int x,y,z,id;}e[N],q[N];bool cmp(P a,P b){return a.z<b.z;} void add(int x,int z){for(;x<N;x+=x&-x)c[x]+=z;} int qu(int x){int sum=0;for(;x;x-=x&-x)sum+=c[x];return sum;} int get(int n,int m){ int i,la,sum=0; for(i=1;i<=min(n,m);i=la+1)sum+=(n/i)*(m/i)*(qu(la=min(n/(n/i),m/(m/i)))-qu(i-1)); return sum&2147483647; } int main(){ for(u[1]=1,i=2;i<N;i++){ if(!f[i])p[++t]=i,u[i]=-1; for(j=1;j<=t&&i*p[j]<N;j++){ f[i*p[j]]=1; if(i%p[j]==0)break; u[i*p[j]]=-u[i]; } } for(i=1;i<N;e[i].x=i++)for(j=i;j<N;j+=i)e[j].z+=i; for(sort(e+1,e+N,cmp),scanf("%d",&T),i=1;i<=T;q[i].id=i++)scanf("%d%d%d",&q[i].x,&q[i].y,&q[i].z); for(sort(q+1,q+T+1,cmp),i=j=1;i<=T;ans[q[i].id]=get(q[i].x,q[i].y),i++) for(;j<N&&e[j].z<=q[i].z;j++)for(k=e[j].x;k<N;k+=e[j].x)add(k,e[j].z*u[k/e[j].x]); for(i=1;i<=T;i++)printf("%d\n",ans[i]); }
BZOJ2154 ans=(i=1~n,j=1~m)Σij/gcd(i,j),令d=gcd(i,j),F(x,y)=(i=1~x,j=1^y,gcd(i,j)=1)Σij
那么可以得到ans=(d=1~min(n,m))Σd^2F(n/d,m/d)/d=ΣdF(n/d)(m/d)
然后对F(x,y)进行反演,F(x,y)=(d=1~min(x,y))Σu(d)Σi(d|i=1~x)Σj(d|j=1~y)=(d=1~min(x,y))Σu(d)d^2Σi(i=1~x/d)Σj(j=1~y/d)=(i=1~min(x,y))Σi^2u(i)(i=1~x,j=1^y)Σij
而很明显2u(i)(i=1~x,j=1^y)Σij=i*(i-1)/2*j*(j-1)/2
这样两个公式都可以n^0.5求出,总时间复杂度O(n)
#include<bits/stdc++.h> #define N 10010000 #define M 20101009 using namespace std;typedef long long LL; int i,j,t,n,m,la,u[N],p[N];LL ans,sum[N];bool f[N]; LL O(LL x){return (x*(x+1)/2)%M;} LL F(int x,int y){ int i,la;LL re=0; for(i=1;i<=min(x,y);i=la+1)re=(re+O(x/i)*O(y/i)%M*(sum[la=min(x/(x/i),y/(y/i))]-sum[i-1]))%M; return re; } int main(){ for(sum[1]=u[1]=1,i=2;i<N;i++){ if(!f[i])p[++t]=i,u[i]=-1; for(j=1;j<=t&&i*p[j]<N;j++){ f[i*p[j]]=1;if(i%p[j]==0)break; u[i*p[j]]=-u[i]; } sum[i]=(sum[i-1]+(LL)i*i*u[i])%M; } for(scanf("%d%d",&n,&m),i=1;i<=min(n,m);i=la+1)ans=(ans+F(n/i,m/i)*(O(la=min(n/(n/i),m/(m/i)))-O(i-1)+M))%M; printf("%lld",ans); }
BZOJ2693 上面一题的多组询问版,考虑继续简化公式
#include<bits/stdc++.h> #define N 10010000 #define M 100000009 using namespace std;typedef long long LL; int i,j,t,n,m,T,la,p[N];LL ans,sum[N],u[N];bool f[N]; LL O(LL x){return (x*(x+1)/2)%M;} int main(){ for(sum[1]=u[1]=1,i=2;i<N;i++){ if(!f[i])p[++t]=i,u[i]=(i-(LL)i*i)%M; for(j=1;j<=t&&i*p[j]<N;j++){ f[i*p[j]]=1; if(i%p[j]==0){ u[i*p[j]]=p[j]*u[i]%M; break; } u[i*p[j]]=u[p[j]]*u[i]%M; } sum[i]=(sum[i-1]+u[i])%M; } for(scanf("%d",&T);T--;printf("%lld\n",(ans+M)%M),ans=0) for(scanf("%d%d",&n,&m),i=1;i<=min(n,m);i=la+1) ans=(ans+O(n/i)*O(m/i)%M*(sum[la=min(n/(n/i),m/(m/i))]-sum[i-1]))%M; }
BZOJ3944 打BC的时候最后时候推出莫比乌斯发现n=1e9没法快速前缀和,记得看到过BZOJ上一道题可以O(n^2/3)求莫比乌斯前缀和,就向Claris问到题号抄程序贴上去A了。比赛完再把这题看了下,大概就是≤n^2/3的预处理,大于n^2/3的递归后MAP快速出解。
#include<bits/stdc++.h> #define M 4000000 using namespace std;typedef long long LL; int n,i,j,x,p,T,pr[M],u[M];bool f[M];LL phi[M];map<int,LL>mp,mo; LL P(int n){ if(n<M)return phi[n];if(mp.count(n))return mp[n]; LL x=1ll*n*(n+1)/2; for(int i=2,la;i<=n;i=la+1)la=n/(n/i),x-=P(n/i)*(la-i+1); return mp[n]=x; } LL U(int n){ if(n<M)return u[n]; if(mo.count(n))return mo[n]; LL x=1,i=2,la; for(i=2,la;i<=n;i=la+1)la=n/(n/i),x-=U(n/i)*(la-i+1); return mo[n]=x; } int main(){ for(phi[1]=u[1]=1,i=2;i<M;i++){ if(!f[i])pr[++p]=i,phi[i]=i-1,u[i]=-1; for(j=1;j<=p&&i*pr[j]<M;j++){ f[i*pr[j]]=1; if(i%pr[j]==0){ u[pr[j]*i]=0; phi[i*pr[j]]=phi[i]*pr[j]; break; } phi[i*pr[j]]=phi[i]*(pr[j]-1); u[pr[j]*i]=-u[i]; } } for(i=2;i<M;i++)u[i]+=u[i-1],phi[i]+=phi[i-1]; for(scanf("%lld",&T);T--;printf("%lld %lld\n",P(x),U(x)))scanf("%d",&x); }
BZOJ4407 ans=Σ(i<=n)Σ(j<=m)gcd(i,j)^k=Σd^kΣ(i<=n/d)Σ(j<=m/d)gcd(i,j)=1
=Σd^kΣu(w)*(n/wd)*(m/wd)=Σ(n/D)(m/D)Σ(d|D)d^k*(u(D/d))
由于快速幂Log消不掉,就暴力求吧。。
#include<bits/stdc++.h> #define M 1000000007 #define N 5000001 using namespace std;typedef long long LL;bool f[N];int T,n,m,k,i,j,t,u[N],p[N];LL ans,v,F[N],h[N]; LL pow(LL a,LL b,LL p){for(v=1;b;a=a*a%p,b>>=1)if(b&1)v=v*a%p;return v;} int main(){ for(scanf("%d%d",&T,&k),i=2;i<N;i++){ h[i]=pow(i,k,M);if(!f[i])p[++t]=i,F[i]=h[i]-1; for(j=1;j<=t&&i*p[j]<N;j++){ f[i*p[j]]=1; if(i%p[j])F[i*p[j]]=(F[i]*(h[p[j]]-1))%M;else{F[i*p[j]]=F[i]*h[p[j]]%M;break;} } } for(F[1]=1,i=2;i<N;i++)F[i]=(F[i]+F[i-1])%M; for(;T--;printf("%lld\n",ans))for(scanf("%d%d",&n,&m),ans=0,i=1;i<=min(n,m);i=j+1)j=min(n/(n/i),m/(m/i)),ans=(ans+(F[j]-F[i-1]+M)*(n/i)%M*(m/i))%M; }
BZOJ3309 枚举d=gcd(i,j),得到
。
如果求出的前缀和就可以O()得到每组答案。
令T=p1a1*p2a2*…*pkak,d=p1b1*p2b2*…*pkbk,0<=ai-bi<=1。
发现如果存在ai!=aj,那么g(T)=0。否则g(T)=(-1)k+1。
BZOJ3512 设。那么在时,有。
两边求和得
对于的情况,找到最大的k|n使得,。
特殊地,。对较小的进行预处理,然后记忆化搜索得到答案。
#include<bits/stdc++.h> #define N 4500000 #define M 1000000007 #define pm make_pair using namespace std;typedef long long LL;LL w,ans;map<pair<LL,LL>, LL>mp; bool f[N];int i,j,t,n,m,u[N],phi[N],p[N],to[N],sm[N]; LL P(LL n,LL m){ if(!m)return 0;if(n==1&&m<N)return sm[m]; if(mp.count(pm(n,m)))return mp[pm(n,m)]; LL i,x,la; if(n==1){ for(x=(m*(m+1)/2)%M,i=2;i<=m;i=la+1)la=m/(m/i),x=(x-P(1,m/i)*(la-i+1)%M+M)%M; return mp[pm(n,m)]=x; } if(u[n]){ for(i=1,x=0;i*i<=n;i++)if(n%i==0)if(x=(x+phi[n/i]*P(i,m/i))%M,i*i!=n)x=(x+phi[i]*P(n/i,m/(n/i)))%M; return mp[pm(n,m)]=x; } return (n/to[n])*P(to[n],m)%M; } int main(){ for(p[0]=phi[1]=u[1]=1,i=2;i<N;i++){ if(!f[i])p[++t]=i,phi[i]=i-1,u[i]=-1; for(j=1;j<=t&&i*p[j]<N;j++){ f[i*p[j]]=1; if(i%p[j]==0){u[p[j]*i]=0;phi[i*p[j]]=phi[i]*p[j];break;} phi[i*p[j]]=phi[i]*(p[j]-1);u[p[j]*i]=-u[i]; } } for(i=1;i<=100000;i++)if(u[i]==0){ for(j=1;j*j<=i;j++)if(i%j==0)if(u[i/j]){to[i]=i/j;break;} if(!to[i])for(j=1;j*j<=i;j++)if(i%j==0)if(u[j])to[i]=j; } for(i=1;i<N;i++)sm[i]=(phi[i]+sm[i-1])%M; for(scanf("%d%d",&n,&m),i=1;i<=n;i++)ans=(ans+P(i,m))%M; printf("%lld\n",(ans+M)%M); }
BZOJ3561 设n<=m,=。
枚举d,然后暴力得到和,时间复杂度O(nlgn)。
#include<bits/stdc++.h> #define W 505050 #define P 1000000007 using namespace std;typedef long long LL;int n,m,i,j,t,d,N,M,p[W],u[W];bool f[W];LL A,w,o,v[W],V[W]; LL pow(LL a,LL b,LL p){for(o=1;b;a=a*a%p,b>>=1)if(b&1)o=o*a%p;return o;} int main(){ scanf("%d%d",&n,&m);if(n>m)swap(n,m); for(u[1]=1,i=2;i<=n;i++){ if(!f[i])p[++t]=i,u[i]=-1; for(j=1;j<=t&&i*p[j]<=n;j++){f[i*p[j]]=1;if(i%p[j]==0)break;u[i*p[j]]=-u[i];} } for(i=1;i<=m;i++)v[i]=1; for(d=1;d<=n;A=(A+w*pow(d,d,P))%P,d++){ for(N=n/d,M=m/d,i=1;i<=M;i++)v[i]=v[i]*i%P,V[i]=(v[i]+V[i-1])%P; for(w=0,i=1;i<=N;i++)if(u[i])w=(w+v[i]*v[i]%P*V[N/i]%P*V[M/i]%P*u[i]+P)%P; }printf("%lld",A); }
待做:3434、3601
2015年6月16日 07:52
莫比乌斯函数只是把nlogn的大力容斥变成了O(n)而已。。
不过如果有修改什么的他非常有用