太阳神拉很喜欢最小公倍数，有一天他想到了一个关于最小公倍数的题目
求满足如下条件的数对(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}}\rfloor}\sum_{j'=1}^{\lfloor{\frac{n}{i'·d}}\rfloor}[gcd(i',j')=1]\]
引入莫比乌斯函数：
\[\sum_{d=1}^{n}\sum_{i'=1}^{\lfloor{\frac{n}{d}}\rfloor}\sum_{j'=1}^{\lfloor{\frac{n}{i'·d}}\rfloor}\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''}}\rfloor\]
把\(d'\)拿出来：
\[\sum_{d'=1}^{\sqrt{n}}\mu(d')[d·i''·j''\leq \frac{n}{d'^2}\text{的组数}]\]
然后就是求满足\(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\n",ans);
return 0;
}
/*
10000000000
210705255
*/原文地址：https://www.cnblogs.com/tqr06/p/11619827.html