[bzoj3994]约数个数和

前排提醒，本文是蒟蒻写的，极大可能性出错，请谨慎食用

好久不写数学题了，随便做一下吧。

题目大意：
给定N,M(<=50000),求：$\sum_{i=1}^{N}\sum_{j=1}^{M}d(ij)$

口胡：
首先我们有一个结论：$\sum_{i=1}^{N}\sum_{j=1}^{M}d(ij)=\sum_{i=1}^{N}\sum_{j=1}^{M}\lfloor\frac{N}{i}\rfloor\lfloor\frac{M}{j}\rfloor[gcd(i,j)==1]$
那么如何证明呢？
首先我们令$F(N,M)=\sum_{i=1}^{N}\sum_{j=1}^{M}d(ij)$
令$G(N,M)=\sum_{i=1}^{N}\sum_{j=1}^{M}\lfloor\frac{N}{i}\rfloor\lfloor\frac{M}{j}\rfloor[gcd(i,j)==1]$
做一个二维差分：

$$F(N,M)-F(N-1,M)-F(N,M-1)+F(N-1,M-1)=d(ij)$$
于是问题转化为了我们要证明： $d(ij)=\sum_{i|N}\sum_{j|M}1*[gcd(i,j)==1]$
我们对于所有素数分类讨论，可以发现一个素数p，假设存在 $n=n'*p^x,m=m'*p^y$那么对于式子左侧我们可以发现这个素数对于答案的贡献就是 $x+y+1$，而式子右侧表示的其实就是 $1<=i<=N,1<=j<=m$并且 $gcd(i,j)==1$的有序数对的个数，那么我们考虑p，右边式子的贡献就是 $(p^x,1),(p^{x-1},1).....(p,1),(1,1),(1,p),(1,p^1).....(1,p^y)$，也只有 $x+y+1$种数对。

于是式子得到了证明：

$$\sum_{i=1}^{N}\sum_{j=1}^{M}d(ij)=\sum_{i=1}^{N}\sum_{j=1}^{M}\lfloor\frac{N}{i}\rfloor\lfloor\frac{M}{j}\rfloor[gcd(i,j)==1]$$
下面就是一个裸的反演啦~
$$\sum_{i=1}^{N}\sum_{j=1}^{M}\lfloor\frac{N}{i}\rfloor\lfloor\frac{M}{j}\rfloor[gcd(i,j)==1]\\ =\sum_{i=1}^{N}\sum_{j=1}^{M}\lfloor\frac{N}{i}\rfloor\lfloor\frac{M}{j}\rfloor\sum_{d|gcd(i,j)}{\mu(d)}\\ =\sum_{d=1}^{N}{\mu(d)}\sum_{i=1}^{\lfloor\frac{N}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{M}{d}\rfloor}{\lfloor\frac{N}{id}\rfloor\lfloor\frac{M}{jd}\rfloor}\\=\sum_{d=1}^{N}{\mu(d)}\sum_{i=1}^{\lfloor\frac{N}{d}\rfloor}{\lfloor\frac{N}{id}\rfloor}\sum_{j=1}^{\lfloor\frac{M}{d}\rfloor}{\lfloor\frac{M}{jd}\rfloor}$$
我们令 $G(N)=\sum_{i=1}^{N}{\lfloor\frac{N}{i}\rfloor}$
于是原式可变形为： $\sum_{d=1}^{N}{\mu(d)}G(\lfloor\frac{N}{d}\rfloor)G(\lfloor\frac{M}{d}\rfloor)$
那么 $G(n)$函数就是群众喜闻乐见的 $d(N)$的前缀和，线性筛约数个数即可。
线性筛复杂度： $O(n)$分块求和复杂度： $O(n^{0.5})$
code：

#include<bits/stdc++.h>
using namespace std;
const int maxn=5e4+10;
typedef long long ll;
ll prime[maxn],mu[maxn];
bool isprime[maxn];
ll c[maxn],d[maxn];
int cnt;
//c[i]表示i的最小因子的次数 
inline void get(){
    mu[1]=c[1]=d[1]=1;
    for(int i=2;i<=maxn;i++){
        if(!isprime[i]){
            prime[++cnt]=i;
            mu[i]=-1;
            c[i]=1;
            d[i]=2;
        }
        for(int j=1;j<=cnt&&i*prime[j]<maxn;j++){
            isprime[i*prime[j]]=true;
            if(i%prime[j]==0){
                d[i*prime[j]]=d[i]*(c[i]+2)/(c[i]+1);
                c[i*prime[j]]=c[i]+1;
                break;
            }
            mu[i*prime[j]]=-mu[i];
            d[i*prime[j]]=d[i]*d[prime[j]];
            c[i*prime[j]]=1;
        }
    }
    for(int i=1;i<maxn;i++)
        mu[i]+=mu[i-1];
    for(int j=1;j<maxn;j++)
        d[j]+=d[j-1];
}
inline int read(){
    int x=0;char ch=getchar();
    while(!isdigit(ch)) ch=getchar();
    while(isdigit(ch)) x=x*10+ch-48,ch=getchar();
    return x;
}
int n,m;
int main(int argc,const char * argv[]){
    get();
    int T=read();
    while(T--){
        n=read();m=read();
        if(n>m) swap(n,m);
        ll ans=0;
        for(int i=1,nex;i<=n;i=nex+1){
            nex=min(n/(n/i),m/(m/i));
            ans+=(mu[nex]-mu[i-1])*d[n/i]*d[m/i];
        }
        printf("%lld\n",ans);
    }
    return 0;
}