Cptraser/ 十月 26, 2018/ 2018.10/ 0 comments

数论分块与莫比乌斯反演(二)

相关链接:数论分块与莫比乌斯反演(一)
题目传送门

d(x)x的约数个数,给定N、M,求\sum_{i=1}^N\sum_{j=1}^M d(ij)
关于d(ij)有以下定理

    \[d(ij)=\sum_{x|i}\sum_{y|j}[gcd(x,y)==1]\]

Tips:证明考虑指数互质,然后感性理解

然后原式就转化成了:

    \[\sum_{i=1}^n\sum_{j=1}^m\sum_{x|i}\sum_{y|j}[gcd(x,y)==1]\]

我们考虑先枚举xy

    \[\sum_{x=1}^n\sum_{y=1}^m\lfloor\frac{n}{x}\rfloor\lfloor\frac{m}{y}\rfloor[gcd(x,y)==1]\]

接着推(日常展开gcd(x,y)==1):

    \[\sum_{x=1}^n\sum_{y=1}^m\lfloor\frac{n}{x}\rfloor\lfloor\frac{m}{y}\rfloor\sum_{d|gcd(x,y)}\mu(d)\]

然后提出d

    \[\sum_{d=1}^{min(n,m)}\mu(d)\sum_{x=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{y=1}^{\lfloor\frac{m}{d}\rfloor}\lfloor\frac{n}{xd}\rfloor\lfloor\frac{m}{yd}\rfloor\]

Tips:x=x'd,y=y'd,其实这里的x,y就是枚举x',y',这也是为什么\lfloor\frac{n}{xd}\rfloor\lfloor\frac{m}{yd}\rfloor的原因,其实它本身的只没有改变,只是枚举的东西改变了。由于x,y提出了d,故范围除以d

再提出每个和式中的常数项:

    \[(\sum_{d=1}^{min(n,m)}\mu(d))(\sum_{x=1}^{\lfloor\frac{n}{d}\rfloor}\lfloor\frac{n}{xd}\rfloor)(\sum_{y=1}^{\lfloor\frac{m}{d}\rfloor}\lfloor\frac{m}{yd}\rfloor)\]

g(x)=\sum_{i=1}^n\lfloor\frac{n}{i}\rfloor,然后原式就变成:

    \[\sum_{d=1}^{min(n,m)}\mu(d)\sum_{x=1}^{\lfloor\frac{n}{d}\rfloor}g(\lfloor\frac{n}{x}\rfloor)\sum_{y=1}^{\lfloor\frac{m}{d}\rfloor}g(\lfloor\frac{m}{y}\rfloor)\]

g O(n\sqrt n)预处理即可。
该式用数论分块求解。

#include <cstdio>
#define V 50005
using namespace std;

int min(int x,int y){return x<y?x:y;}

long long Ans=0,F[V];
int N,M,T,Cnt,K;
int Prime[V],vis[V],Mul[V],Sum[V];
void init(void){
    Mul[1]=Sum[1]=1;
    for(int i=2;i<=V-5;i++){
        if(!vis[i])Prime[++Cnt]=i,Mul[i]=-1;
        for(int j=1;j<=Cnt;j++){
            if(Prime[j]*i>V)break;
            vis[Prime[j]*i]=1;
            if(!(i%Prime[j]))break;
            Mul[i*Prime[j]]=-Mul[i];
        }
        Sum[i]=Sum[i-1]+Mul[i];
    }
    for(int i=1;i<=V;i++){
        for(int L=1,R;L<=i;L=R+1){
            R=i/(i/L);
            F[i]+=1ll*(R-L+1)*(i/L);
        }
    }
}


int main()
{
    for(init(),scanf("%d",&T);T;T--){
        scanf("%d%d",&N,&M),K=min(N,M),Ans=0;
        for(int L=1,R;L<=K;L=R+1){
            R=min(N/(N/L),M/(M/L));
            Ans+=1ll*(Sum[R]-Sum[L-1])*F[N/L]*F[M/L];
        }
        printf("%lld\n",Ans);
    }
}
Share this Post

Leave a Comment

电子邮件地址不会被公开。 必填项已用*标注

You may use these HTML tags and attributes: <a href="" title=""> <abbr title=""> <acronym title=""> <b> <blockquote cite=""> <cite> <code> <del datetime=""> <em> <i> <q cite=""> <s> <strike> <strong>
*
*