[bzoj3994] [SDOI2015] 约数个数和

Description

设 \(d(x)\) 为 \(x\) 的约数个数，给定 \(N\) 、\(M\)，求 \(\sum\limits_{i=1}^N \sum\limits_{j=1}^M d(ij)\)

Input

输入文件包含多组测试数据。
第一行，一个整数 \(T\) ，表示测试数据的组数。
接下来的 \(T\) 行，每行两个整数 \(N\) 、\(M\) 。

Output

\(T\) 行，每行一个整数，表示你所求的答案。

Sample Input

2

7 4

5 6

Sample Output

110

121

HINT

\(1<=N, M<=50000\)
\(1<=T<=50000\)

---

想法

首先有一个结论
$ d(i,j)=\sum\limits_{x|i} \sum\limits_{y|j} [gcd(x,y)==1] $
原谅我不会证 \(qwq\)

接下来，喜闻乐见的推式子

\[\begin{equation*} \begin{aligned} &\sum\limits_{i=1}^n \sum\limits_{j=1}^m d(ij) \\ =&\sum\limits_{i=1}^n \sum\limits_{j=1}^m \sum\limits_{x|i} \sum\limits_{y|j} [gcd(x,y)==1] \\ =&\sum\limits_{i=1}^n \sum\limits_{j=1}^m [gcd(x,y)==1] \times \lfloor \frac{n}{x} \rfloor \times \lfloor \frac{m}{y} \rfloor \\ \end{aligned} \end{equation*} \]

设 \(f(i)=\sum\limits_{i=1}^n \sum\limits_{j=1}^m [gcd(x,y)==1] \times \lfloor \frac{n}{x} \rfloor \times \lfloor \frac{m}{y} \rfloor\)
应用莫比乌斯反演

\[\begin{equation*} \begin{aligned} F(i)&=f(i)+f(2i)+... \\ &=\sum\limits_{x=1}^{\lfloor \frac{n}{i} \rfloor} \sum\limits_{y=1}^{\lfloor \frac{m}{i} \rfloor} \lfloor \frac{n}{ix} \rfloor \lfloor \frac{m}{iy} \rfloor \\ &=\sum\limits_{x=1}^{\lfloor \frac{n}{i} \rfloor} \lfloor \frac{n}{ix} \rfloor \sum\limits_{y=1}^{\lfloor \frac{m}{i} \rfloor} \lfloor \frac{m}{iy} \rfloor \end{aligned} \end{equation*} \]

不妨设 \(s(x)=\sum\limits_{i=1}^x \lfloor \frac{x}{i} \rfloor\)
它可以预处理后 \(O(1)\) 调用
则 \(f(i)=\sum\limits_{i|d} \mu(\frac{d}{i}) s(\frac{n}{d}) s(\frac{m}{d})\)
那么

\[\begin{equation*} \begin{aligned} ans&=f(1) \\ &=\sum\limits_{i=1}^{min(n,m)} \mu(i) s(\frac{n}{i})s(\frac{m}{i}) \end{aligned} \end{equation*} \]

然后就好啦。

---

代码

#include<cstdio>
#include<iostream>
#include<algorithm>
 
using namespace std;
 
const int N = 50005;
typedef long long ll;
 
ll s[N];
int mu[N],p[N],prime[N],pnum,sum[N];
void getp(){
    mu[1]=sum[1]=1;
    for(int i=2;i<N;i++) p[i]=1;
    for(int i=2;i<N;i++){
        if(p[i]) prime[pnum++]=i,mu[i]=-1;
        for(int j=0;j<pnum && (ll)prime[j]*i<N;j++){
            p[i*prime[j]]=0;
            if(i%prime[j]==0){
                mu[i*prime[j]]=0;
                break;
            }
            mu[i*prime[j]]=-mu[i];
        }
        sum[i]=sum[i-1]+mu[i];
    }
}
void gets(){
    for(int x=1;x<N;x++){
        for(int l=1,r;l<=x;l=r+1){
            r=x/(x/l);
            s[x]+=1ll*(x/l)*(r-l+1);
        }
    }
}
 
int main()
{
    int T,n,m;
    getp(); gets();
     
    scanf("%d",&T);
    while(T--){
        scanf("%d%d",&n,&m);
        ll ans=0;
        for(int l=1,r;l<=n && l<=m;l=r+1){
            r=min(n/(n/l),m/(m/l));
            ans+=s[n/l]*s[m/l]*(sum[r]-sum[l-1]);
        }
        printf("%lld\n",ans);
    }
     
    return 0;
}