莫比乌斯反演

形式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 于神之怒加强版

click here！ 

#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 数表

click here！

#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

click here!

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