莫比乌斯反演
形式1
已经有函数F(n)=∑ f(d)，可以导出 f(n)= ∑ μ(d)F(n/d)
d|n d|n
形式2
已经有F(n)= ∑ f(d),可以导出f(n)= ∑ μ(d/n)F(d)
n|d n|d
bzoj2301 Problem b
题目大意：求gcd（x,y）=k（a<=x<=b,c<=y<=d）的数对个数。
思路：F（i）表示i|gcd（x,y）的数对个数，f(i)表示gcd（x，y）的数对个数，我们要求的就是f（k）。F（i）相对好求一些，F（i）=⌊n/i⌋⌊m/i⌋，
f(i)=∑（i|a）μ(a/i)F(a)=∑（i|a）μ(a|i)⌊n/a⌋⌊m/a⌋ ，然后我们只需要枚举⌊n/a⌋⌊m/a⌋的值就可以了（这个值据说是2（根n+根m））。我们穷举a，
找到所有⌊n/a⌋⌊m/a⌋相等的的位置（连续的），然后乘上预处理出来的mu函数的和就可以了。这里有一个小技巧，1...x的数列中，与n/i相等的值的最
后一位是（n/(n/i)），这样我们找到n、m右边界较小的就是一段相等的区间了。
同时，对于求一个区间内gcd=k的，我们可以通过将区间/k后找到互质的数对的个数就是答案了。
再同时，对于这道题中的答案要容斥原理一下，每次都处理成[1...x][1...y]，然后答案就是work（b，d）-work（a，d）-work（c，b）+work（a，c）了。
这题中的技巧和化简思路都十分重要，这种穷举除数的思路据（某大神）说十分常用。
#include<iostream> #include<cstdio> #include<cstring> #include<algorithm> using namespace std; int mu[50010]={0},prime[50010]={0},k; bool flag[50010]={false}; void prework(int n) { int i,j; mu[1]=1; for (i=2;i<=n;++i) { if (!flag[i]) { prime[++prime[0]]=i;mu[i]=-1; } for (j=1;j<=prime[0]&&i*prime[j]<=n;++j) { flag[i*prime[j]]=true; if (i%prime[j]==0) { mu[i*prime[j]]=0;break; } else mu[i*prime[j]]=-mu[i]; } } for (i=1;i<=n;++i) mu[i]+=mu[i-1]; } int work(int x,int y) { int i,last=0,ans=0; x/=k;y/=k; if (x>y) swap(x,y); for (i=1;i<=x;i=last+1) { last=min(x/(x/i),y/(y/i)); ans+=(x/i)*(y/i)*(mu[last]-mu[i-1]); } return ans; } int main() { int n,a,b,c,d,i,ans; prework(50000); scanf("%d",&n); for (i=1;i<=n;++i) { scanf("%d%d%d%d%d",&a,&b,&c,&d,&k); ans=work(b,d)-work(a-1,d)-work(c-1,b)+work(a-1,c-1); printf("%d\n",ans); } }
bzoj2818 GCD
题目大意：求1～n中gcd(x,y)为素数的对数。
思路：可以筛出素数，然后求出n/prime[i]内的gcd(x,y)=1的数对就可以了。
#include<iostream> #include<cstdio> #include<cstring> #include<algorithm> #define maxnode 10000005 #define LL long long using namespace std; int prime[maxnode]={0},mu[maxnode]={0}; bool flag[maxnode]={false}; void shai(int n) { int i,j;mu[1]=1; for (i=2;i<=n;++i) { if (!flag[i]) {prime[++prime[0]]=i;mu[i]=-1;} for (j=1;j<=prime[0]&&prime[j]*i<=n;++j) { flag[prime[j]*i]=true; if (i%prime[j]==0){mu[i*prime[j]]=0;break;} else mu[i*prime[j]]=-mu[i]; } } for (i=1;i<=n;++i) mu[i]+=mu[i-1]; } LL work(int x) { int i,last;LL ans=0; for (i=1;i<=x;i=last+1) { last=x/(x/i); ans+=(LL)(x/i)*(LL)(x/i)*(LL)(mu[last]-mu[i-1]); } return ans; } int main() { int n,i,j;LL ans=0; scanf("%d",&n);shai(n); for (i=1;i<=prime[0];++i) ans+=work(n/prime[i]); printf("%I64d\n",ans); }
当然也可以求一个phi的前缀和sum，然后1+2*sum[n/prime[i]]，也是答案。
bzoj3853GCD Array
题目大意：在一个数列上进行操作。1）读入p，d，v，对所有gcd（i,p）=d的ai+v；2）求$\sum_{i=1}^xa_i$.
思路：听了题解才明白是反演。我们对于一次修改操作可以看作$a_i+=[gcd(i,p)=d]v$，对上式进行化简可以得到
[gcd(i,p)=d]v=[gcd (\frac{i}{d},\frac{p}{d})=1]v
=v\sum_{e|\frac{i}{d},e| \frac{p}{d} } \mu{e}
= v\sum_{ed|i,e|\frac{p}{d}}\mu{e}
=\sum_{ed|i , e| \frac{p}{d}} \mu{e}v
我们维护一个数组$f_e$，令$a_i=\sum_{e|i}f_e$，修改的时候我们每次穷举$\frac{p}{d}$的约数x（也就是e），x*d的位置加上x*v，用树状数组维护；查询的时候ANS=\sum_{i=1}^xa_i
=\sum_{i=1}^x\sum_{e|i}f(e)
=\sum_{i=1}^x\left \lfloor\frac{x}{i}\right \rfloor f(i)
#include<iostream> #include<cstdio> #include<cstring> #include<algorithm> #include<vector> #define maxnode 200005 using namespace std; long long cc[maxnode]={0}; int n,prime[maxnode]={0},mu[maxnode]={0}; vector<int> yue[maxnode]; bool flag[maxnode]={false}; void pre() { int i,j; mu[1]=1; for (i=2;i<maxnode;++i) { if (!flag[i]) { prime[++prime[0]]=i;mu[i]=-1; } for (j=1;prime[j]*i<maxnode&&j<=prime[0];++j) { flag[prime[j]*i]=true; if (i%prime[j]) mu[i*prime[j]]=-mu[i]; else { mu[i*prime[j]]=0;break; } } } } void getyue() { int i,j; for (i=1;i<=maxnode;++i) yue[i].clear(); for (i=1;i<maxnode;++i) for (j=i;j<maxnode;j+=i) yue[j].push_back(i); } int lowbit(int x){return x&(-x);} long long ask(int x,int y) { long long sum=0; --x; for (;y;y-=lowbit(y)) sum+=cc[y]; for (;x;x-=lowbit(x)) sum-=cc[x]; return sum; } void change(int x,int v) { for (;x<=n;x+=lowbit(x)) cc[x]+=(long long)v; } int main() { int m,i,j,p,v,d,op,ci=0; long long ans; pre();getyue(); while(scanf("%d%d",&n,&m)) { if (n==0&&m==0) break; memset(cc,0,sizeof(cc));++ci; printf("Case #%d:\n",ci); for (i=1;i<=m;++i) { scanf("%d",&op); if (op==1) { scanf("%d%d%d",&p,&d,&v); if (p%d) continue; p=p/d; for (j=0;j<yue[p].size();++j) change(yue[p][j]*d,mu[yue[p][j]]*v); } else { scanf("%d",&p);ans=0; for (j=1;j<=p;j=d+1) { d=(p/(p/j)); ans+=(long long)(p/j)*ask(j,d); } printf("%lld\n",ans); } } } }
bzoj2005 noi2010能量采集
题目大意：n×m的矩阵中，有n×m个点位于格点上，对于每个点，采集的能量是这个点和（0，0）连线上点的个数（不包含端点）×2+1，求采集所有点的能量和。
思路：对于一个点(x,y)，采集的能量是2×（gcd（x,y）-1)+1，求这个的sigma。可以用反演，穷举1～min（n，m）作为gcd的值，求出个数后乘i加起来乘2-1就是答案了。（反演的过程同第一题）。
其实并不需要这么麻烦，我们在反演中F（d）表示d|gcd（x，y）的个数，F(d)=(n/d)(m/d)，那么我们用这个数减去gcd=2d、3d……id的个数就可以了，用f（d）表示gcd=d的个数，就是用F（d）-f（2d）-f（3d）-……f（id），我们倒着穷举gcd，就可以完成了。不需要反演和求mu，速度还++。。。
#include<iostream> #include<cstdio> #include<cstring> #include<algorithm> #define maxnode 100005 #define up 100000 using namespace std; int mu[maxnode]={0},prime[maxnode]={0}; long long sum[maxnode]={0}; bool flag[maxnode]={0}; void getmu() { int i,j; mu[1]=1; for (i=2;i<=up;++i) { if (!flag[i]) { prime[++prime[0]]=i;mu[i]=-1; } for (j=1;j<=prime[0]&&i*prime[j]<=up;++j) { flag[i*prime[j]]=true; if (i%prime[j]==0){mu[i*prime[j]]=0;break;} else mu[i*prime[j]]=-mu[i]; } } for (i=1;i<=up;++i) sum[i]=sum[i-1]+(long long)mu[i]; } long long ask(int n,int m,int k) { int i,j,last=0; long long ans=0; n/=k;m/=k;if (n>m) swap(n,m); for (i=1;i<=n;i=last+1) { last=min(n/(n/i),m/(m/i)); ans+=(long long)(n/i)*(long long)(m/i)*(sum[last]-sum[i-1]); } return ans; } int main() { int n,m,i,j; long long ans=0; scanf("%d%d",&n,&m); if (n>m) swap(n,m);getmu(); for (i=1;i<=n;++i) ans+=i*ask(n,m,i); ans=ans*2-(long long)n*(long long)m; printf("%I64d\n",ans); }
在暑假跟着faebdc犇学习莫比乌斯反演，发现只需要记住两个式子，剩下的就是熟练转换了。
sigma phi(d)=n sigma mu(d)=e(n) ===>e(i)= (i==1 ? 1 : 0)
d|n d|n
bzoj2154 Crash的数字表格
题目大意：求sigma（i=1～n）sigma（j=1～m） lcm（i，j）
思路；原式=sigma sigma i*j/gcd(i,j)
i=1~n j=1~m
=sigma sigma k1*k2*gcd(i,j) (设i=k1gcd（i，j），j=k2gcd（i，j）,gcd(k1,k2)=1)
i=1~n j=1~m
= sigma d sigma sigma i*j*e(gcd(i,j))
d=1~min(n,m) i=1~n/d j=1~m/d
= sigma d sigma sigma i*j sigma mu(k)
d=1~min(n,m) i=1~n/d j=1~m/d k|i&k|j
= sigma d sigma k^2 mu(k) sum(n/d/k,m/d/k) sum(x,y)=sigma sigma i*j=(i*(i+1)/2)*(j*(j+1)/2)
d=1~min(n,m) k=1~min(n/d,m/d) i=1~x j=1~y
根号套根号的求出两个sigma，O（n）复杂度。（注意强转，i,j互质 <==> e(lca(i,j))=1）
#include<iostream> #include<cstring> #include<algorithm> #include<cstdio> #define maxnode 10000005 #define p 20101009 #define LL long long using namespace std; int prime[maxnode]={0},mu[maxnode]={0}; LL sum[maxnode]={0}; bool flag[maxnode]={false}; void shai(int n) { int i,j; mu[1]=1; for (i=2;i<=n;++i) { if (!flag[i]) { prime[++prime[0]]=i;mu[i]=-1; } for (j=1;j<=prime[0]&&i*prime[j]<=n;++j) { flag[i*prime[j]]=true; if (i%prime[j]==0) { mu[i*prime[j]]=0;break; } else mu[i*prime[j]]=-mu[i]; } } for (i=1;i<=n;++i) sum[i]=(sum[i-1]+(((LL)mu[i]*(((LL)i*(LL)i)%p))%p+p)%p)%p; } LL getsum(LL *a,int l,int r) { --l;return ((a[r]-a[l])%p+p)%p; } LL get(int n,int m) { int i,last=0;LL ans=0; if (n>m) swap(n,m); for (i=1;i<=n;i=last+1) { last=min(n/(n/i),m/(m/i)); LL xx=((LL)(n/i)*(LL)(n/i+1)/2%p)*((LL)(m/i)*(LL)(m/i+1)/2%p)%p; ans=(ans+getsum(sum,i,last)*xx%p)%p; } return (ans%p+p)%p; } LL work(int n,int m) { int i,last=0;LL ans=0; if (n>m) swap(n,m); for (i=1;i<=n;i=last+1) { last=min(n/(n/i),m/(m/i)); ans=(ans+((LL)(last+i)*(LL)(last-i+1)/2%p)*get(n/i,m/i))%p; } return (ans%p+p)%p; } int main() { int n,m,i,j; scanf("%d%d",&n,&m);shai(max(n,m)); printf("%lld",work(n,m)); }
bzoj2671 Calc
题目大意：求1~n中，(a+b)|ab的个数。
思路：因为之前做过一道求ab=k(a+b)的问题，所以总想从这方面下手，但是思考了很久都没有进展。后来改变了一下思路，其实可以将式子都除去一个gcd(a,b)，得到(x+y)|xygcd(a,b)，因为gcd(x+y,x)=1,gcd(x+y,y)=1（orz sunshine）,所以(x+y)|gcd(a,b)，那么会给答案贡献n/(y*(x+y))（一）（这个式子是因为设gcd(a,b)=t(x+y)，会有b=t(x+y)y<=n，所以有这么多个给答案。注意分母不为零！！！），其中gcd(x,y)=1,x<y，这里的y的上界是根n(因为gcd(a,b)*b<=n)，所以我们枚举一下y，x可以n^(1/4)枚举，对于每一段（一）相等的段中，求出与y互质的数的个数，这里用到了容斥原理和莫比乌斯函数，一段区间[l,r]中与一个数n互质的数的个数是sigma(k|n)mu[k]*(r/k-(l-1)/k)，所以就可以算出答案了。去网上找了一下这种做法的复杂度，暴力出奇迹吧。
#include<iostream> #include<cstdio> #include<cstring> #include<algorithm> #include<cmath> #define maxnode 100005 #define LL long long using namespace std; int prime[maxnode]={0},mu[maxnode]={0},yue[maxnode*10]={0}; bool flag[maxnode]={false}; void shai(int n) { int i,j;mu[1]=1; for (i=2;i<n;++i){ if (!flag[i]){prime[++prime[0]]=i;mu[i]=-1;} for (j=1;j<=prime[0]&&prime[j]*i<n;++j){ flag[prime[j]*i]=true; if (i%prime[j])mu[i*prime[j]]=-mu[i]; else{mu[i*prime[j]]=0;break;} } } } void fen(LL n) { LL i;yue[0]=0; for (i=1;i*i<=n;++i) if (n%i==0){ yue[++yue[0]]=i; if (i*i!=n) yue[++yue[0]]=n/i; } } LL work(LL n) { LL a,b,l,r,i,j,up,last,ans=0;up=sqrt(n); for (b=1;b<=up;++b){ fen(b); for (a=1;a<b&&b*(a+b)<=n;a=last+1){ last=min(n/(n/(b*(a+b)))/b-b,b-1); for (i=1;i<=yue[0];++i) ans+=n/(b*(a+b))*mu[yue[i]]*(last/yue[i]-(a-1)/yue[i]); } }return ans; } int main() { int n;scanf("%d",&n);shai(sqrt(n)+1); printf("%I64d\n",work((LL)n)); }
bzoj3994 约数个数和
题目大意：设d(x)表示x的约数个数，求sigma(i=1~n)sigma(j=1~m)d(ij)
思路：d(ij)=sigma(i|n)sigma(j|m)e(gcd(i,j))（考虑如果i、j不互质的时候可以通过质因数的变换使得互质，这样就导致多次计算了。）对这个式子反演。
sigma(i=1~n)sigma(j=1~m)d(ij)=sigma(i=1~n)sigma(j=1~m)(n/i)(m/j)e(gcd(i,j))
=sigma(i=1~n)sigma(j=1~m)(n/i)(m/j)sigma(d|i&&d|j)mu(d)
=sigma(d=1~min(n,m))mu(d)gi(n/d)gi(m/d) (gi(x)=sigma(i=1~x)(x/i)，这一步可以设i=k1gcd(),j=k2gcd() 考虑)
#include<iostream> #include<cstdio> #include<cstring> #include<algorithm> #define maxm 50005 #define LL long long using namespace std; int prime[maxm]={0},mu[maxm]={0}; LL gi[maxm]={0}; bool flag[maxm]={0}; void shai(int n){ int i,j;mu[1]=1; for (i=2;i<=n;++i){ if (!flag[i]){prime[++prime[0]]=i;mu[i]=-1;} for (j=1;j<=prime[0]&&prime[j]*i<=n;++j){ flag[prime[j]*i]=true; if (i%prime[j]) mu[i*prime[j]]=-mu[i]; else{mu[i*prime[j]]=0;break;} } }for (i=2;i<=n;++i) mu[i]+=mu[i-1]; } void pre(int n){ int i,j,la; for (i=1;i<=n;++i) for (j=1;j<=i;j=la+1){ la=i/(i/j);gi[i]+=(LL)(i/j)*(LL)(la-j+1); } } LL calc(int n,int m){ int i,j,la;LL ans=0; if (n>m) swap(n,m); for (i=1;i<=n;i=la+1){ la=min(n/(n/i),m/(m/i)); ans+=(LL)(mu[la]-mu[i-1])*gi[n/i]*gi[m/i]; }return ans; } int main(){ int t,n,m; shai(maxm-1);pre(maxm-1); scanf("%d",&t); while(t--){ scanf("%d%d",&n,&m); printf("%I64d\n",calc(n,m)); } }
bzoj3561 DZY Loves Math VI
题目大意：求sigma(i=1~n)sigma(j=1~m) lcm(i,j)^gcd(i,j)。
思路：sigma(i=1~n)sigma(j=1~m) lcm(i,j)^gcd(i,j)=sigma(i=1~n)sigma(j=1~m)[k1*k2*gcd(i,j)]^gcd(i,j)*e(gcd(k1,k2))
=sigma(d=1~min(n,m))d^d sigma(a=1~min(n/d,m/d)) a^(2d)*mu(a)g(n/d/a,m/d/a,d)
设g(nn,mm,d)=sigma(i=1~nn)i^d * sigma(j=1~mm) j^d，这个数组可以用vector预处理，是nlnn的。
做的时候暴力循环就可以了，复杂度是nlnn的。
#include<iostream> #include<cstdio> #include<cstring> #include<algorithm> #include<vector> #define N 500005 #define p 1000000007LL #define LL long long using namespace std; vector<LL> sum[N]; int prime[N]={0},mu[N]; LL mi[N]; bool flag[N]={false}; void shai(){ int i,j;mu[1]=1; for (i=2;i<N;++i){ if (!flag[i]){ prime[++prime[0]]=i;mu[i]=-1; }for (j=1;j<=prime[0]&&i*prime[j]<N;++j){ flag[i*prime[j]]=true; if (i%prime[j]) mu[i*prime[j]]=-mu[i]; else{mu[i*prime[j]]=0;break;} } } } LL getg(int x,int y,int d){return sum[d][x]*sum[d][y]%p;} LL getmi(LL x,int y){ LL a=1LL; for (;y;y>>=1){ if (y&1) a=a*x%p; x=x*x%p; }return a;} int main(){ int i,j,n,m,up;LL ans=0LL,ci; scanf("%d%d",&n,&m); if (n>m) swap(n,m); shai(); for (i=1;i<=m;++i){mi[i]=1LL;sum[i].push_back(0LL);} for (i=1;i<=m;++i) for (j=1;j*i<=m;++j){ mi[j]=mi[j]*(LL)j%p; sum[i].push_back((sum[i][j-1]+mi[j])%p); } for (i=1;i<=m;++i) mi[i]=1LL; for (i=1;i<=n;++i){ up=min(n/i,m/i); for (ci=0LL,j=1;j<=up;++j){ mi[j]=mi[j]*(LL)j%p; ci=(ci+mi[j]*mi[j]%p*mu[j]*getg(n/i/j,m/i/j,i))%p; }ci=(ci%p+p)%p; ans=(ans+getmi((LL)i,i)*ci%p)%p; }printf("%I64d\n",ans); }
bzoj4407 于神之怒加强版
#include<iostream> #include<cstring> #include<cstdio> #include<algorithm> #define N 5000005 #define LL long long #define p 1000000007 using namespace std; int k,prime[N]={0}; LL gi[N]={0},ci[N]; bool flag[N]={0}; LL mi(LL x,int y){ LL a=1LL; for(;y;y>>=1){ if (y&1) a=a*x%p; x=x*x%p; }return a;} void shai(){ int i,j,v;gi[1]=1; for (i=2;i<N;++i){ if (!flag[i]){ prime[++prime[0]]=i; ci[i]=mi((LL)i,k); gi[i]=((ci[i]-1LL)%p+p)%p; }for (j=1;j<=prime[0]&&i*prime[j]<N;++j){ flag[v=i*prime[j]]=true; if (i%prime[j]) gi[v]=(gi[i]*(ci[prime[j]]-1)%p+p)%p; else{ gi[v]=gi[i]*ci[prime[j]]%p; break;} } }for (i=2;i<N;++i) gi[i]=(gi[i]+gi[i-1])%p;} int main(){ int t,n,m,i,j,up,last;LL ans; scanf("%d%d",&t,&k);shai(); while(t--){ scanf("%d%d",&n,&m);ans=0LL; for (up=min(n,m),i=1;i<=up;i=last+1){ last=min(n/(n/i),m/(m/i)); ans=(ans+(LL)(n/i)*(LL)(m/i)%p*((gi[last]-gi[i-1])%p+p)%p)%p; }printf("%I64d\n",ans); } }
bzoj3529 数表
#include<iostream> #include<cstdio> #include<cstring> #include<algorithm> #define N 100005 #define LL long long #define p 2147483648LL using namespace std; struct use{int n,m,a,po;LL ans;}ai[N]={0}; struct uu{LL v;int po;}fi[N]={0}; LL ci[N]={0LL}; int prime[N]={0},flag[N]={0},cnt[N]={0},mu[N]={0}; int in(){ int x=0;char ch=getchar(); while(ch<'0'||ch>'9') ch=getchar(); while(ch>='0'&&ch<='9'){ x=x*10+ch-'0';ch=getchar(); }return x;} LL mi(LL x,int y){ if (!y) return 1LL; if (y==1) return x; LL mm=mi(x,y/2); if (y%2) return mm*mm*x; else return mm*mm;} void shai(){ int i,j,x,y;LL xx;mu[1]=1;fi[1].v=1LL; for (i=1;i<N;++i) fi[i].po=i; for (i=2;i<N;++i){ if (!flag[i]){ mu[i]=-1;fi[i].v=(LL)i+1; prime[++prime[0]]=i;cnt[i]=1; }for (j=1;j<=prime[0]&&(x=i*prime[j])<N;++j){ flag[x]=1; if (i%prime[j]){ fi[x].v=fi[i].v*(LL)(prime[j]+1); mu[x]=-mu[i];cnt[x]=1; }else{ cnt[x]=cnt[i]+1;xx=mi(prime[j],cnt[x]); fi[x].v=fi[i].v*(xx*(LL)prime[j]-1LL)/(xx-1LL); mu[x]=0;} } } } int lowbit(int x){return x&(-x);} int cmp(const uu&x,const uu&y){return x.v<y.v;} int cmp1(const use&x,const use&y){return x.a<y.a;} int cmp2(const use&x,const use&y){return x.po<y.po;} void tch(int x,LL y){for (;x<N;x+=lowbit(x)) ci[x]=(ci[x]+y)%p;} void insert(int x){ int i,j,y;y=fi[x].po; for (i=(N-1)/y;i;--i) if (mu[i]!=0) tch(i*y,(LL)mu[i]*fi[x].v);} LL ask(int x){ LL sum=0LL; for (;x;x-=lowbit(x)) sum=(sum+ci[x])%p; return sum;} LL sig(int l,int r){return ((ask(r)-ask(l-1))%p+p)%p;} int main(){ int t,i,j,x,y;t=in();shai(); for (i=1;i<=t;++i){ai[i].n=in();ai[i].m=in();ai[i].a=in();ai[i].po=i;} sort(ai+1,ai+t+1,cmp1); sort(fi+1,fi+N,cmp); for (j=1,i=1;i<=t;++i){ while(fi[j].v<=ai[i].a&&j<N) insert(j++); if (ai[i].n>ai[i].m) swap(ai[i].n,ai[i].m); for (x=1;x<=ai[i].n;x=y+1){ y=min(ai[i].n/(ai[i].n/x),ai[i].m/(ai[i].m/x)); ai[i].ans=(ai[i].ans+(LL)(ai[i].m/x)*(LL)(ai[i].n/x)%p*sig(x,y)%p)%p; } }sort(ai+1,ai+t+1,cmp2); for (i=1;i<=t;++i) printf("%I64d\n",ai[i].ans); }
bzoj3309 DZY Loves Math
#include<iostream> #include<cstdio> #include<cstring> #include<algorithm> #define N 10000000 #define LL long long using namespace std; int prime[N+1]={0},fi[N+1]={0}; bool flag[N+1]={false}; void dfs(int i,int k,int ci){ if (i>prime[0]) return; if ((LL)ci*(LL)prime[i]>N) return; LL v=(LL)(ci*prime[i]); dfs(i+1,-k,ci*prime[i]); for (LL vv=v;vv<=N;vv*=v) fi[vv]=k; dfs(i+1,k,ci);} void shai(){ int i,j,v;fi[0]=fi[1]=0; for (i=2;i<=N;++i){ if (!flag[i]) prime[++prime[0]]=i; for (j=1;j<=prime[0]&&i*prime[j]<=N;++j){ flag[v=i*prime[j]]=true; if (i%prime[j]==0) break; } }dfs(1,1,1); for (i=2;i<=N;++i) fi[i]+=fi[i-1]; } int main(){ int t,a,b,i,la;LL ans; shai();scanf("%d",&t); while(t--){ scanf("%d%d",&a,&b); if (a>b) swap(a,b); for (ans=0LL,i=1;i<=a;i=la+1){ la=min(a/(a/i),b/(b/i)); ans+=(LL)(a/i)*(LL)(b/i)*(LL)(fi[la]-fi[i-1]); }printf("%I64d\n",ans); } }
bzoj3202 项链
click here!（！！！）
#include<iostream> #include<cstdio> #include<cstring> #include<algorithm> #define N 10000005 #define LL long long #define pp 1000000007LL using namespace std; int prime[N]={0},mu[N],a,bi[N],az; LL inv,ai[N],n,ans,smu[N],p; bool flag[N]={false}; void shai(){ int i,j,v;mu[1]=1;smu[1]=1LL; for (i=2;i<N;++i){ if (!flag[i]){prime[++prime[0]]=i;mu[i]=-1;} for (j=1;j<=prime[0]&&(v=i*prime[j])<N;++j){ flag[v]=true; if (i%prime[j]) mu[v]=-mu[i]; else{mu[v]=0;break;} }smu[i]=smu[i-1]+(LL)mu[i]; } } LL mul(LL x,LL y){ LL aa=0LL;x=(x%p+p)%p;y=(y%p+p)%p; if (y>x) swap(x,y); for (;y;y>>=1LL){ if (y&1LL) aa=(aa+x)%p; x=(x+x)%p; }return aa;} LL mi(LL x,LL y){ LL aa=1LL;x=(x%p+p)%p; for (;y;y>>=1LL){ if (y&1LL) aa=mul(aa,x)%p; x=mul(x,x)%p; }return aa;} LL getc(LL x){return mul(mul(mul(x,x+1LL),(x+2LL)),inv);} LL calc(){ LL an=0LL;int i,la; for (i=1;i<=a;i=la+1){ la=a/(a/i); an=(an+mul(smu[la]-smu[i-1],getc((LL)(a/i))))%p; }return (an%p+p)%p;} LL getf(LL nn,LL x){ return ((mi(x-1LL,nn)+(nn%2LL ? -1LL : 1LL)*(x-1LL)%p)%p+p)%p; } void dfs(int i,LL x,LL y,LL z){ if (i==az+1){ ans=((ans+mul(getf(n/x,z),y))%p+p)%p; return; }int j;LL cc=1LL;dfs(i+1,x,y,z); for (j=1;j<=bi[i];++j){ cc*=ai[i]; dfs(i+1,x*cc,y*(cc/ai[i]*(ai[i]-1LL)),z); } } int main(){ LL x;int t,i; shai();scanf("%d",&t); while(t--){ scanf("%I64d%d",&n,&a);az=0; for (x=n,i=1;i<=prime[0]&&x>1LL;++i) if (x%(LL)prime[i]==0LL){ ai[++az]=(LL)prime[i];bi[az]=0; while(x%(LL)prime[i]==0LL){++bi[az];x/=(LL)prime[i];} } if (x>1){ai[++az]=x;bi[az]=1;} if (n%pp==0LL) p=pp*pp; else p=pp; inv=mi(6LL,pp*(pp-1LL)-1LL); x=calc();ans=0LL;dfs(1,1LL,1LL,x); if (n%pp==0LL){p=pp;ans=ans/p%p*mi(n/p,p-2LL)%p;} else ans=mul(ans,mi(n,p-2LL)); printf("%I64d\n",ans); } }
