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;
}

__EOF__

本文作者：acwing_zyy
本文链接：https://www.cnblogs.com/zyyun/p/16423901.html
关于博主：评论和私信会在第一时间回复。或者直接私信我。
版权声明：本博客所有文章除特别声明外，均采用 BY-NC-SA 许可协议。转载请注明出处！
声援博主：如果您觉得文章对您有帮助，可以点击文章右下角【推荐】一下。您的鼓励是博主的最大动力！