【莫比乌斯反演】[HAOI2011]Problem b

题目链接

P2522 [HAOI2011]Problem b

题目大意

求：

$$\sum _{i=x}^n\sum _{j=y}^m[\gcd (i,j)=k]$$

题解

首先这是一个二维前缀和的形式，我们改写一下：

$$\sum _{i=1}^n\sum _{j=1}^m[\gcd (i,j)=k]-\sum _{i=1}^{x-1}\sum _{j=1}^m[\gcd (i,j)=k]-\sum _{i=1}^n\sum _{j=1}^{y-1}[\gcd (i,j)=k]+\sum _{i=1}^{x-1}\sum _{j=1}^{y-1}[\gcd (i,j)=k]$$

定义 $s(n,m)=\sum _{i=1}^n\sum _{j=1}^m[\gcd (i,j)=k]$，直接考虑求解 $s(n,m)$。

枚举每一种 $\gcd (i,j)$，令 $i=i^{\prime }k,j=j^{\prime }k$，得到：

$$\sum _{i=1}^{\lfloor \frac{n}{k}\rfloor }\sum _{j=1}^{\lfloor \frac{m}{k}\rfloor }[\gcd (i,j)=1]$$

利用莫比乌斯反演 $[n=1]=\sum _{d\mid n}\mu (d)$，得到：

$$\sum _{i=1}^{\lfloor \frac{n}{k}\rfloor }\sum _{j=1}^{\lfloor \frac{m}{k}\rfloor }\sum _{d\mid \gcd (i,j)}\mu (d)$$

让 $d$ 枚举 $\gcd (i,j)$ 的每一个因数，得：

$$\sum _{d=1}^n\sum _{i=1}^{\lfloor \frac{n}{k}\rfloor }\sum _{j=1}^{\lfloor \frac{m}{k}\rfloor }\mu (d)\times [d\mid i\text{and}d\mid j]$$

将所有式子各自归为：

$$\sum _{d=1}^n\mu (d)\sum _{i=1}^{\lfloor \frac{n}{k}\rfloor }[d\mid i]\sum _{j=1}^{\lfloor \frac{m}{k}\rfloor }[d\mid j]$$

这个时候继续照常处理，提出 $d$：

$$\sum _{d=1}^n\mu (d)\sum _{i=1}^{\lfloor \frac{n}{kd}\rfloor }[1\mid i]\sum _{j=1}^{\lfloor \frac{m}{kd}\rfloor }[1\mid j]$$

所以式子等于：

$$\sum _{d=1}^n\mu (d)\left(\lfloor \frac{n}{kd}\rfloor \lfloor \frac{m}{kd}\rfloor \right)$$

因此整除分块即可，时间复杂度 $\mathcal{O}(\sqrt{n})$。

代码

//你 cnt=1 了吗？
#include<bits/stdc++.h>
#define mid ((l+r)>>1)
#define fir first
#define sec second
#define lowbit(i) (i&(-i))
using namespace std;
const int N=5e4+5;
const int inf=1e18;
struct edge{int to,nxt,l;};
inline int read(){
    char op=getchar();
    int w=0,s=1;
    while(op<'0'||op>'9'){
        if(op=='-') s=-1;
        op=getchar();
    }
    while(op>='0'&&op<='9'){
        w=(w<<1)+(w<<3)+op-'0';
        op=getchar();
    }
    return w*s;
}
//数论部分
const double pi=acos(-1);
const int mod=1e18-11;
int Mul(int a,int b){return (a%mod*b%mod)%mod;}
int Add(int a,int b){return (a+b)%mod;}
int Dec(int a,int b){return (a-b+mod)%mod;}
int Pow(int a,int k){
    int ans=1;
    while(k){
        if(k&1) ans=Mul(ans,a);
        a=Mul(a,a);
        k>>=1;
    }
    return ans;
}
int gcd(int x,int y){return y==0?x:gcd(y,x%y);}
int lcm(int x,int y){return x/gcd(x,y)*y;}
int inv(int x){return Pow(x,mod-2);}
void exgcd(int a,int b,int &x,int &y){
    if(b==0){
        x=1,y=0;
        return;
    }
    exgcd(b,a%b,x,y);
    int t=x;
    x=y;
    y=t-a/b*y;
}
int n,m,mu[N],vis[N],p[N],sum[N],cnt=0,k;
void mobius(){
    mu[1]=1;
    for(register int i=2;i<=N-5;i++){
        if(!vis[i]){
            mu[i]=-1;
            p[++cnt]=i;
        }
        for(register int j=1;j<=cnt&&p[j]*i<=N-5;j++){
            vis[i*p[j]]=1;
            if(i%p[j]==0){
                mu[i*p[j]]=0;
                break;
            }
            mu[i*p[j]]=-mu[i];
        }
    }
    for(register int i=1;i<=N-5;i++) sum[i]=sum[i-1]+mu[i];
}
long long Sqrt(int n,int m){
    if(!n||!m) return 0;
    int l=1,r;
    long long ans=0;
    for(l=1;l<=min(n,m);l=r+1){
        r=min(n/(n/l),m/(m/l));
        ans+=(long long)(sum[r]-sum[l-1])*(n/(l*k))*(m/(l*k));
    }
    return ans;
}
signed main(){
    mobius();
    int T=read();
    while(T--){
        int x=read(),n=read(),y=read(),m=read();
        k=read();
        printf("%lld\n",Sqrt(n,m)-Sqrt(x-1,m)-Sqrt(n,y-1)+Sqrt(x-1,y-1));
    }
}