[20191003机房测试] 太阳神

太阳神拉很喜欢最小公倍数,有一天他想到了一个关于最小公倍数的题目
求满足如下条件的数对(a,b)对数:
a,b 均为正整数且 a,b<=n 而lcm(a,b)>n
其中的 lcm 当然表示最小公倍数
答案对 1,000,000,007取模

数据范围是1e10的,打表找了半天规律发现没用……
那就莫比乌斯反演呗

题目求:

[sum_{i=1}^{n}sum_{j=1}^{n}[lcm(i,j)> n] ]

但是大于的太多了,那我们就反过来求,最后用总的来减

也就是求:

[sum_{i=1}^{n}sum_{j=1}^{n}[lcm(i,j)leq n] ]

先把(lcm)拆开

[sum_{i=1}^{n}sum_{j=1}^{n}[dfrac{i·j}{gcd(i,j)}leq n] ]

看到(gcd),枚举 (i,j) 的约数 (d)

[sum_{d=1}^{n}sum_{i=1}^{n}sum_{j=1}^{n}[dfrac{i·j}{d}leq n] ]

拉出去反演:

[sum_{d=1}^{n}sum_{i'=1}^{lfloor{frac{n}{d}} floor}sum_{j'=1}^{lfloor{frac{n}{i'·d}} floor}[gcd(i',j')=1] ]

引入莫比乌斯函数:

[sum_{d=1}^{n}sum_{i'=1}^{lfloor{frac{n}{d}} floor}sum_{j'=1}^{lfloor{frac{n}{i'·d}} floor}mu(gcd(i',j')) ]

再枚举 (i',j') 的约数 (d'),继续反演

[sum_{d=1}^{n}sum_{d'=1}^{sqrt{frac{n}{d}}}mu(d')sum_{i''=1}sum_{j''=1}lfloor{frac{n}{d·d'^2i''}} floor ]

(d')拿出来:

[sum_{d'=1}^{sqrt{n}}mu(d')[d·i''·j''leq frac{n}{d'^2} ext{的组数}] ]

然后就是求满足(d·i''·j''leq frac{n}{d'^2})的三元组(d,i'',j'')的个数

这个可以(Theta(n^{frac{2}{3}}))算出来

就可以了

代码:

#include<bits/stdc++.h>
#define ll long long
#define N 100005
#define M 100000
#define mod 1000000007
using namespace std;

ll n,ans,siz,sum,maxn,x;
ll mu[N],prime[N],primenum;
bool isprime[N];

template<class T>inline void read(T &res)
{
	char c;T flag=1;
	while((c=getchar())<'0'||c>'9')if(c=='-')flag=-1;res=c-'0';
	while((c=getchar())>='0'&&c<='9')res=res*10+c-'0';res*=flag;
}

void init()
{
	mu[1]=1;
	for(register int i=2;i<=M;++i)
	{
		if(!isprime[i])
		{
			prime[++primenum]=i;
			mu[i]=-1;
		}
		for(register int j=1;j<=primenum;++j)
		{
			if(i*prime[j]>M) break;
			isprime[i*prime[j]]=1;
			if(!(i%prime[j]))
			{
				mu[i*prime[j]]=0;
				break;
			}
			mu[i*prime[j]]=-mu[i];
		}
	}
}

int main()
{
	freopen("ra.in","r",stdin);
	freopen("ra.out","w",stdout);
	read(n);
	init();
	siz=(ll)(sqrt(n)+0.5);
	for(register ll i=1;i<=siz;++i)
	{
		sum=0;
		maxn=n/i/i;
		for(register ll j=1;j*j*j<=maxn;++j)
		{ 
			for(register ll k=j;k*k<=maxn/j;++k)
			{
				x=(maxn/j/k-k+1);
				if(j==k) sum=(sum+1+(x-1)*3)%mod;
				else sum=(sum+3+(x-1)*6)%mod;
			}
		}
		sum%=mod;
		ans=(ans+mu[i]*sum)%mod;
	}
	ans=((n%mod)*(n%mod)-ans)%mod;
	while(ans<0) ans+=mod;
	printf("%lld
",ans);
	return 0;
}
/*
10000000000
210705255
*/
原文地址:https://www.cnblogs.com/tqr06/p/11619827.html