51nod 1222 最小公倍数计数

题意

对于区间，我们求出([1,b])的答案减去([1,a-1])的即可。

求([1,n])的答案：
(ans=sumlimits_{i=1}^{n}sumlimits_{j=1}^n[frac{i*j}{gcd(i,j)}leqslant n])
(=sumlimits_{d=1}^nsumlimits_{i=1}^{frac{n}{d}}sumlimits_{j=1}^{frac{n}{d}}[i*jleqslant frac{n}{d}][gcd(i,j)=1])
(=sumlimits_{d=1}^nsumlimits_{k=1}^{frac{n}{d}}mu(k)sumlimits_{i=1}^{frac{n}{dk}}sumlimits_{j=1}^{frac{n}{dk}}[i*k*j*kleqslant frac{n}{d}])
(=sumlimits_{d=1}^nsumlimits_{k=1}^{frac{n}{d}}mu(k)sumlimits_{i=1}^{frac{n}{dk}}sumlimits_{j=1}^{frac{n}{dk}}[i*j*dleqslant frac{n}{k^2}])
(=sumlimits_{k=1}^nmu(k)sumlimits_{d=1}^{frac{n}{k}}sumlimits_{i=1}^{frac{n}{dk}}sumlimits_{j=1}^{frac{n}{dk}}[i*j*dleqslant frac{n}{k^2}])
(=sumlimits_{k=1}^{sqrt{n}}mu(k)sumlimits_{d=1}^{frac{n}{k^2}}sumlimits_{i=1}^{frac{n}{dk^2}}sumlimits_{j=1}^{frac{n}{dk^2}}[i*j*dleqslant frac{n}{k^2}])
于是我们枚举(k)，算出满足([d*i*jleqslant frac{n}{k^2}])的三元组((d,i,j))有多少。

我们发现((d,i,j))的上界似乎没什么用，我们只要找到一个三元组((a,b,c))满足(a*b*cleqslant frac{n}{k^2})，就一定能满足上面的关系。

于是我们就是要求三元组((a,b,c))满足(a*b*cleqslant frac{n}{k^2})的个数。

我们分别统计三个都不相同，有两个相同和三个都相同的方案数。

我们先求各不相同的，枚举小的那个，设为(x)，(x^3leqslant frac{n}{k^2})。之后枚举第二大的，设为(y)，(x*y*yleqslant frac{n}{k^2})，最大的那个的范围就能确定。中途再算一下((x,x,y),(x,y,y),(x,x,x))这种，统计一下就好了。

code:

#include<bits/stdc++.h>
using namespace std;
#define int long long
const int maxn=1e6; 
int a,b;
int mu[maxn];
bool vis[maxn];
vector<int>prime;
inline int read()
{
	char c=getchar();int res=0,f=1;
	while(c<'0'||c>'9'){if(c=='-')f=-1;c=getchar();}
	while(c>='0'&&c<='9')res=res*10+c-'0',c=getchar();
	return res*f;
}
inline void prework(int n)
{
	vis[1]=mu[1]=1;
	for(int i=2;i<=n;i++)
	{
		if(!vis[i])prime.push_back(i),mu[i]=-1;
		for(unsigned int j=0;j<prime.size()&&i*prime[j]<=n;j++)
		{
			vis[i*prime[j]]=1;
			if(i%prime[j]==0)break;
			mu[i*prime[j]]=mu[i]*mu[prime[j]];
		}
	}
}
inline int calc(int n)
{
	int res=0;
	for(int k=1;k*k<=n;k++)
	{
		if(!mu[k])continue;
		int tmp=0,up=n/(k*k);
		for(int i=1;i*i*i<=up;i++)
		{
			for(int j=i+1;j*j*i<=up;j++)tmp+=(up/(i*j)-j)*6+3;//三个不同+(i,j,j)
			tmp+=(up/(i*i)-i)*3;//(i,i,j)
			tmp++;//(i,i,i)
		}
		res+=mu[k]*tmp;
	}
	return (res+n)/2;
}
signed main()
{
	prework(1e6);
	a=read(),b=read();
	printf("%lld",calc(b)-calc(a-1));
	return 0;
}