51nod-1222-最小公倍数计数

题意

给到 \(a,b\) ，求

\[\sum _{i=a}^b\sum _x\sum _y[x\le y][\text{lcm}(x,y)=i] \]

即最小公倍数在 \([a,b]\) 中的有序数对个数。\(a,b\le 10^{11}\) 。

分析

转化成求 \(\sum _{x}\sum _{y}[\text{lcm}(x,y)\le n]\) ，最后加上 \(x=y\) 的情况除以2即可得到有序数对的个数。减一减即可得到答案。

\[\begin{aligned} \sum _{x=1}^n\sum _{y=1}^n[\frac{xy}{\gcd(x,y)}\le n]&=\sum _{d=1}^n\sum _{x=1}^{\lfloor\frac{n}{d}\rfloor}\sum _{y=1}^{\lfloor\frac{n}{d}\rfloor}[xyd\le n][\gcd(x,y)=1] \\ &=\sum _{d=1}^n\sum _{x=1}^{\lfloor\frac{n}{d}\rfloor}\sum _{y=1}^{\lfloor\frac{n}{d}\rfloor}[xyd\le n]\sum _{e|x,e|y}\mu (e) \\ &=\sum _{d=1}^n\sum _{e=1}^{\lfloor\frac{n}{d}\rfloor}\mu (e)\sum _{x=1}^{\lfloor\frac{n}{de}\rfloor}\sum _{y=1}^{\lfloor\frac{n}{de}\rfloor}[xyde^2\le n] \\ &=\sum _{e=1}^n\mu (e)\sum _{d=1}^{\lfloor\frac{n}{e}\rfloor}\sum _{x=1}^{\lfloor\frac{n}{de}\rfloor}\sum _{y=1}^{\lfloor\frac{n}{de}\rfloor}[xyd\le \lfloor\frac{n}{e^2}\rfloor] \\ \end{aligned} \]

注意到后面是个 \(\lfloor\frac{n}{e^2}\rfloor\) ，所以 \(e\) 只需要枚举到 \(\sqrt n\) 。

\[ans=\sum _{e=1}^\sqrt n\mu (e)\sum _{d=1}^{\lfloor\frac{n}{e}\rfloor}\sum _{x=1}^{\lfloor\frac{n}{de}\rfloor}\sum _{y=1}^{\lfloor\frac{n}{de}\rfloor}[xyd\le \lfloor\frac{n}{e^2}\rfloor] \\ \]

讨论求和上界。若 \(d>\lfloor\frac{n}{e}\rfloor\) ，那么一定不满足 \(d\le \lfloor\frac{n}{e^2}\rfloor\) 。如果 \(x\ge \lfloor\frac{n}{de}\rfloor\) ，那么一定不可能有 \(xd\ge \lfloor\frac{n}{e^2}\rfloor\) ，因此右边的三个求和上界都是无效的，可以被后面的条件直接限制。

于是现在问题就变成了求

\[f(n)=\sum _x\sum _y\sum _z[xyz\le n] \]

注意到三个求和项是相同的，地位相等可以轮换的，不如强行给它们定序。设 \(x\le y\le z\) ，计算完后乘上 6 即可扩展到所有排列。再讨论一下 \(x=y\le z,x\le y=z,x=y=z\) 的情况，容斥一下即可算出答案。

现在考虑如何算第一个 \(x\le y\le z\) 的情况。由于我们有了序，所以 \(x\le \sqrt [3]n\) ，只需要枚举这些 \(x\) 。枚举完 \(x\) ，剩下的有 \(y\le \sqrt \frac{n}{x}\) ，最后 \(z\) 的取值个数可以直接 \(\lfloor\frac{n}{ij}\rfloor-j+1\) 计算出来。剩下的三种情况都可以 \(O(\sqrt [3]n)\) 或 \(O(1)\) 解决，因此直接拿第一种情况来分析复杂度。

对于 \(f(n)\) ，它的复杂度为

\[\int _0^\sqrt [3]n\sqrt \frac{n}{x} dx=O(n^\frac{2}{3}) \]

对于外层，复杂度为

\[\int _1^\sqrt nf(\frac{n}{x^2})dx=O(n^{\frac{2}{3}}) \]

所以是一个奇妙的暴力。

代码

#include<bits/stdc++.h>
using namespace std;
typedef long long giant;
const int maxn=1e6+1;
bool np[maxn];
int p[maxn],ps=0;
giant mu[maxn];
giant les(giant n) {
	giant fir=0,sec=0,thr=0;
	for (giant i=1;i*i*i<=n;++i) {
		giant ed=(giant)sqrt((long double)n/i);
		for (giant j=i;j<=ed;++j) fir+=n/(i*j)-j+1;
		sec+=n/(i*i)-i+1;
		sec+=(giant)sqrt((long double)n/i)-i+1;
		thr=i;
	}
	return fir*6-sec*3+thr;
}
giant calc(giant n) {
	giant ret=0;
	for (giant e=1,tmp;(tmp=e*e)<=n;++e) 
		ret+=les(n/tmp)*mu[e];
	return (ret+=n)>>=1;	
}
int main() {
#ifndef ONLINE_JUDGE
	freopen("test.in","r",stdin);
#endif
	mu[1]=1;
	for (int i=2;i<maxn;++i) {
		if (!np[i]) p[++ps]=i,mu[i]=-1;
		for (int j=1,tmp;j<=ps && (tmp=i*p[j])<maxn;++j) {
			np[tmp]=true;
			if (i%p[j]==0) break;
			mu[tmp]=-mu[i];
		}
	}
	giant l,r;
	cin>>l>>r;
	giant ans=calc(r)-calc(l-1);
	cout<<ans<<endl;
	return 0;
}