1358. 约数个数和
题目链接
1358. 约数个数和
设 \(d(x)\) 为 \(x\) 的约数个数，给定 \(N,M\)，求
\(\sum_{i=1}^N\sum_{j=1}^Md(ij)\)
输入格式
输入多组测试数据。
第一行，一个整数 \(T\)，表示测试数据的组数。
接下来的 \(T\) 行，每行两个整数 \(N、M\)。
输出格式
\(T\) 行，每行一个整数，表示你所求的答案。
数据范围
\(1 \le N,M,T \le 50000\)
输入样例：
2
7 4
5 6
输出样例：
110
121
解题思路
莫比乌斯反演
结论：\(d(i j)=\sum_{x \mid i} \sum_{y \mid j}[\operatorname{gcd}(x, y)=1]\)
证明：对 \(i,j\) 分别分解质因数：\(i=p_{1}^{\alpha_{1}} \cdots p_{k}^{\alpha_{k}}, j=p_{1}^{\beta_{1}} \cdots p_{k}^{\beta_{k}}, i \times j=p_{1}^{\alpha_{1}+\beta_{1}} \cdots p_{k}^{\alpha_{k}+\beta_{k}}\)，对于 \(p_1\) 来说，要使 \(x \mid i,y \mid j\) 且 \(gcd(x,y) =1\)，则共有 \(\alpha_{1}+\beta_{1}+1\) 种选择，由于 \(p_i\) 之间相互独立，每种 \(p_i\) 有 \(\alpha_{i}+\beta_{i}+1\) 种选择，从这些选择中各选择一个即可组成满足要求的数，共 \(\left(\alpha_{1}+\beta_{1}+1\right) \cdots\left(\alpha_{k}+\beta_{k}+1\right)\) 种方案，该式正好为约数个数
设 \(f(n)=\sum_{i=1}^{N} \sum_{j=1}^{M} \sum_{x \mid i} \sum_{y \mid j}[\operatorname{gcd}(x, y)=n]\)，则由莫比乌斯反演第二条定理，得 \(F(n)=\sum_{i=1}^{N} \sum_{j=1}^{M} \sum_{x \mid i} \sum_{y \mid j}[n \mid \operatorname{gcd}(x, y)]\)，考虑化简 \(F(n)\)，将内层两重循环换到外面，\(x,y\) 范围分别为 \([1,N],[1,M]\)，此时要求 \(x \mid i,y \mid j\)，而 \([n \mid \operatorname{gcd}(x, y)]\) 只与 \(x,y\) 有关，此时范围内 \(x\) 的倍数有 \(N/x\) 个，\(y\) 的倍数有 \(M/y\) 个，\(F(n)=\sum_{x=1}^{N} \sum_{y=1}^{M}(N/x)\times (M/y)\times [n \mid \operatorname{gcd}(x, y)]\)，\(x,y\) 都为 \(n\) 的倍数，设 \(x'=x/n,y'=y/n\)，\(F(n)=\sum_{x'=1}^{N/n} \sum_{y'=1}^{M/n}(N/(x'n))\times (M/(y'n))\)，令 \(N=N/n,M=M/n,x=x',y=y'\)，得 \(F(n)=\sum_{x=1}^{N} \sum_{y=1}^{M}(N/x)\times (M/y)=\sum_{x=1}^{N}(N/x)\times \sum_{y=1}^{M}(M/y)\)，设 \(H(x)=\sum_{i=1}^{x}(x/i)\)，则 未变换之前，\(F(n)=H(N/n)\times H(M/n)\)，而求解 \(H(x)\) 可整除分块 \(O(n\sqrt{n})\) 预处理出来，由莫比乌斯反演，\(f(n)=\sum_{n \mid d} \mu\left(\frac{d}{n}\right) F(d)=\sum_{n \mid d} \mu\left(\frac{d}{n}\right)\times H(N/d)\times H(M/d)\)，而题目要求求解 \(f(1)\)，则 \(f(1)=\sum_{ d} \mu\left(d\right)\times H(N/d)\times H(M/d)\)，其中 \(d\) 为正整数，同时也可整除分块 \(O(\sqrt{n})\) 求解该式
- 时间复杂度：\(O(n\sqrt{n})\)
代码
// Problem: 约数个数和
// Contest: AcWing
// URL: https://www.acwing.com/problem/content/1360/
// Memory Limit: 64 MB
// Time Limit: 2000 ms
//
// Powered by CP Editor (https://cpeditor.org)
// %%%Skyqwq
#include <bits/stdc++.h>
//#define int long long
#define help {cin.tie(NULL); cout.tie(NULL);}
#define pb push_back
#define fi first
#define se second
#define mkp make_pair
using namespace std;
typedef long long LL;
typedef pair<int, int> PII;
typedef pair<LL, LL> PLL;
template <typename T> bool chkMax(T &x, T y) { return (y > x) ? x = y, 1 : 0; }
template <typename T> bool chkMin(T &x, T y) { return (y < x) ? x = y, 1 : 0; }
template <typename T> void inline read(T &x) {
int f = 1; x = 0; char s = getchar();
while (s < '0' || s > '9') { if (s == '-') f = -1; s = getchar(); }
while (s <= '9' && s >= '0') x = x * 10 + (s ^ 48), s = getchar();
x *= f;
}
const int N=5e4+5;
int t,n,m,prime[N],cnt,v[N],u[N],s[N],h[N];
int g(int a,int b)
{
return a/(a/b);
}
void init(int n)
{
u[1]=1;
for(int i=2;i<=n;i++)
{
if(v[i]==0)
{
u[i]=-1;
prime[++cnt]=v[i]=i;
}
for(int j=1;j<=cnt&&i*prime[j]<=n;j++)
{
if(v[i]<prime[j])break;
if(i%prime[j]==0)u[i*prime[j]]=0;
else
u[i*prime[j]]=-u[i];
v[i*prime[j]]=prime[j];
}
}
for(int i=1;i<=n;i++)s[i]=s[i-1]+u[i];
for(int i=1;i<=n;i++)
{
for(int l=1,r;l<=i;l=r+1)
{
r=g(i,l);
h[i]+=(r-l+1)*(i/l);
}
}
}
int main()
{
init(N-1);
cin>>t;
while(t--)
{
cin>>n>>m;
LL res=0;
int sz=min(n,m);
for(int l=1,r;l<=sz;l=r+1)
{
r=min(g(n,l),g(m,l));
res+=1ll*(s[r]-s[l-1])*h[n/l]*h[m/l];
}
cout<<res<<'\n';
}
return 0;
}
