【莫比乌斯反演】[SDOI2015]约数个数和

原题链接

P3327 [SDOI2015]约数个数和

题目大意

求：

$$\sum _{i=1}^n\sum _{j=1}^{n}d(i\times j)$$

其中 $d(i\times j)$ 表示 $i\times j$ 的因数个数。

题解

首先证明公式：

$$d(ij)=\sum _{x\mid i}\sum _{y\mid j}[\gcd (x,y)=1]$$

证：

我们先考虑 $i=p^{k_1},j=p^{k_2}$，其中 $p$ 是质数，现在我们要统计 $ij$ 有多少种不同的因子，那么我们先考虑 $x=1$，继续枚举 $y$，发现在 $j$ 中可以一直选到 $p^{k_2}$ 而做到不重不漏，其中得到的这一堆值分别是 $1,p,p^2,\dots ,p^{k_2}$。

接下来我不能再继续选取上述的因子，考虑怎么选出 $p^{k_2+1}$ 这个因子，显然我只需要让 $x=1$ 变成 $x=p$，此时我们发现对于 $x=p$，我们已经不能再继续选取 $y=p^{0,\dots ,k_2}$，因为已经被 $x=1$ 选过了，所以 $x=p$ 的贡献是 $1$。同理此时我们已经选出了 $p^{k_2+i}$ 这个因子，那我们只有令 $x=p^{i+1},y=p^{k_2}$ 才能得到 $p^{k_2+i+1}$ 这个因子。

所以简而言之，在这种情况下，我们总共能选取的因子个数是 $k_2+k_1-1$。此时我们还可以发现，当我们考虑 $x\ne 1$ 时，由于我们永远选取了 $y=p^{k_2}$，所以我们每个 $x$ 只能贡献一次，这一次我不一定是 $y=p^{k_2}$，如果我们假定 $y=1$，那么 $x\ne 1$ 的贡献仍然是一次，相当于我们令这个 $y=p^{k_2}$ 的贡献转移给了 $y=1$。这个时候我们发现我们所有的贡献都来自于 $x,y$ 中至少取了一个 $1$。也就是 $\gcd (x,y)=1$ 的全部情况。

接下来推广到 $i=\prod p_i^{k_i},j=\prod p_j^{k_j}$。考虑我们在 $i$ 中选出了 $x=1$，此时 $y$ 可以取到任何值，且最多只能取到 $j$。此时如果我们想要取到 $p_{x}j$，就必须依赖 $i$ 中的 $p_x$，而且这个 $p_x$ 不能与任何 $y<j$ 的 $y$ 结合，因为已经被 $x=1$ 考虑过了。所以每一个 $x\ne 1$ 都只会贡献一次 $y=j$，接下来我们不让 $y=j$ 产生贡献，转移到 $y=1$ 与 $x\ne 1$ 产生贡献。因此发现 $x,y$ 中总是至少选取了一个 $1$，故 $\sum _{x\mid i}\sum _{y\mid j}[\gcd (x,y)=1]=d(ij)$。

举个例子方便理解，令 $i=2^3,j=2^5$。我们先考虑 $x=1$ 的情况，此时 $y$ 可以是 $1,2,4,8,16,32$。$ij$ 中肯定有 $64$ 这个因子，只靠 $x=1$ 是枚举不出来的，所以我们将 $x=1$ 改为 $x=2$，这个时候 $y=1,2,4,8,16$ 的情况都会与 $x=1$ 的情况重合，所以我们不能考虑这些 $y$，因此只有 $y=32$ 的时候我们取到了我们以前没有取到的数 $64$。因此对于 $x\ne 1$，都只有 $y=32$ 产生了贡献，这不方便计算，产生了一次贡献可以视作是 $y=1$ 的贡献，尽管乘起来不是原来没有的因子，但反正我产生了一次贡献，因此我至少会在 $x,y$ 中取一个 $1$。对于推广的情况同理。

带入公式得到：

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

转化 $[n=1]$ 的式子：

$$\sum _{i=1}^n\sum _{j=1}^m\sum _{x\mid i}\sum _{y\mid j}\sum _{d\mid \gcd (x,y)}\mu (d)$$

考虑先枚举 $x,y$：

$$\sum _{x=1}^n\sum _{y=1}^m\sum _{i=1}^n\sum _{j=1}^m[x\mid i\text{and}y\mid j]\sum _{d\mid \gcd (x,y)}\mu (d)$$

令 $i=i^{\prime }x,j=j^{\prime }y$，得：

$$\sum _{x=1}^n\sum _{y=1}^m\sum _{i=1}^{\lfloor \frac{n}{x}\rfloor }\sum _{j=1}^{\lfloor \frac{m}{y}\rfloor }[x\mid ix\text{and}y\mid jy]\sum _{d\mid \gcd (x,y)}\mu (d)$$

中括号内式子等于 $1$，我们直接省去：

$$\sum _{x=1}^n\sum _{y=1}^m\sum _{d\mid \gcd (x,y)}\mu (d)\left(\lfloor \frac{n}{x}\rfloor \lfloor \frac{m}{y}\rfloor \right)$$

继续拆掉 $d\mid \gcd (x,y)$：

$$\sum _{d=1}^n\sum _{x=1}^n\sum _{y=1}^m[d\mid x\text{and}d\mid y]\mu (d)\left(\lfloor \frac{n}{x}\rfloor \lfloor \frac{m}{y}\rfloor \right)$$

令 $x=x^{\prime }d,y=y^{\prime }d$，得到：

$$\sum _{d=1}^n\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$$

接下来各回各家：

$$\sum _{d=1}^n\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$$

令 $f(n)$ 表示 $\sum _{i=1}^n\lfloor \frac{n}{i}\rfloor$，得到：

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

内层的 $f$ 函数可以整除分块解决，但为了不让每次询问都进行一次 $\mathcal{O}(\sqrt{n})$ 的运算，我们提前算好这个前缀和。剩下的每次查询也是整除分块，提前预处理 $\mu (d)$ 的前缀和即可，总时间复杂度 $\mathcal{O}(T\sqrt{n})$。

代码

//P3327
//你 cnt=1 了吗？
#include<bits/stdc++.h>
#define int long long
#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=1e9+7;
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 mu[N],vis[N],p[N],cnt=0,ans[N],sum[N];
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&&i*p[j]<=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];
    for(register int i=1;i<=N-5;i++){//预处理整除分块
        for(register int l=1,r;l<=i;l=r+1){
            r=i/(i/l);
            ans[i]+=(i/l)*(r-l+1);
        }
    }
}
int Sqrt(int n,int m){
    int tot=0;
    for(register int l=1,r;l<=min(n,m);l=r+1){
        r=min(n/(n/l),m/(m/l));
        tot+=(sum[r]-sum[l-1])*ans[n/l]*ans[m/l];
    }
    return tot;
}
signed main(){
    int T=read();
    mobius();
    while(T--){
        int n=read(),m=read();
        printf("%lld\n",Sqrt(n,m));
    }
}