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

题意

给到 (a,b) ，求

[sum _{i=a}^bsum _xsum _y[xle y][ ext{lcm}(x,y)=i] ]

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

分析

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

[egin{aligned} sum _{x=1}^nsum _{y=1}^n[frac{xy}{gcd(x,y)}le n]&=sum _{d=1}^nsum _{x=1}^{lfloorfrac{n}{d} floor}sum _{y=1}^{lfloorfrac{n}{d} floor}[xydle n][gcd(x,y)=1] \ &=sum _{d=1}^nsum _{x=1}^{lfloorfrac{n}{d} floor}sum _{y=1}^{lfloorfrac{n}{d} floor}[xydle n]sum _{e|x,e|y}mu (e) \ &=sum _{d=1}^nsum _{e=1}^{lfloorfrac{n}{d} floor}mu (e)sum _{x=1}^{lfloorfrac{n}{de} floor}sum _{y=1}^{lfloorfrac{n}{de} floor}[xyde^2le n] \ &=sum _{e=1}^nmu (e)sum _{d=1}^{lfloorfrac{n}{e} floor}sum _{x=1}^{lfloorfrac{n}{de} floor}sum _{y=1}^{lfloorfrac{n}{de} floor}[xydle lfloorfrac{n}{e^2} floor] \ end{aligned} ]

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

[ans=sum _{e=1}^sqrt nmu (e)sum _{d=1}^{lfloorfrac{n}{e} floor}sum _{x=1}^{lfloorfrac{n}{de} floor}sum _{y=1}^{lfloorfrac{n}{de} floor}[xydle lfloorfrac{n}{e^2} floor] \ ]

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

于是现在问题就变成了求

[f(n)=sum _xsum _ysum _z[xyzle n] ]

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

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

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

[int _0^sqrt [3]nsqrt 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;
}