题目描述
太阳神拉很喜欢最小公倍数,有一天他想到了一个关于最小公倍数的题目。
求满足如下条件的数对(a,b)对数:a,b均为正整数且a,b<=n而lcm(a,b)>n。其中的lcm当然表示最小公倍数。答案对1,000,000,007取模
推式子
我们先做个补集转化。
然后就是统计这样的三元组(a,b,d)个数了
1、a*b*d<=n
2、(a,b)=1
a和b定了,d的取值就有
nab
ans=∑na=1∑nb=1[(a,b)=1]∗nab
我们尝试莫比乌斯反演一下,也就是
设
f[k]=∑na=1∑nb=1[(a,b)=k]∗nab
g[k]=∑nki=1f[k∗i]
反演一下变成
f[k]=∑ni=1μ(i)∗g[k∗i]
g是很容易算的,所以
f[k]=∑ni=1μ(i)∗∑nk∗ia′=1∑nk∗ib′=1(nk2)a′b′
现在我们需要f[1],代入k=1得
ans=∑nd=1μ(d)∗∑nd2a′=1∑nd2b′=1(nd2)a′b′
观察后面,可以改变d的枚举上界
ans=∑n√d=1μ(d)∗∑nd2a′=1∑nd2b′=1(nd2)a′b′
后面那个东西只与
nd2
有关,设其为
f(nd2)
算f
设
f(x)=∑xi=1∑xj=1xij
可以看出可以根据
xi
的取值分块!每一块答案显然相同,块数为
x√
假如i确定了,那么答案为
∑xij=1(xi)j
后面那个东西只与
xi
有关,设其为
g(xi)
算g
设
g(x)=∑xi=1xi
我们应该怎么算g呢?
有这样一个奇妙的方法。
首先,先知道一个结论。
g(x)=∑xi=1d(i)
d函数表示约数个数。
为什么呢?你考虑
xi
,那么说明
1∗i,2∗i,3∗i……xi∗i
都在x内,而且i是上述提到的
xi
个数的约数。
因此可以表示为约数个数前缀和。
那么也就是说g可以线筛?
但显然不能全都线筛啊?
dalao说,以
n2/3
为界,记为N。我们需要筛出g(1~N),然后对于x<=N便直接得解。
如果x>N?观察g的算式,我们显然可以按照取值分块,得到一个根号算法。那我们就这样做最快咯。不过这样会T,加个记忆化即可。
怎么记忆化呢?dalao说,x可以映射到n/x,这样一定不会重复。这样就可以开小记忆化数组。
总复杂度?dalao说,是
O(n2/3)
#include<cstdio>
#include<algorithm>
#include<cmath>
#define fo(i,a,b) for(i=a;i<=b;i++)
using namespace std;
typedef long long ll;
const int N=3000000,mo=1000000007;
int mu[N+10],sum[N+10],mem[N+10],d[N+10],pri[N+10];
bool bz[N+10];
int i,j,k,l,t,m,top,ans;
ll n;
void prepare(){
int i,j;
mu[1]=1;
fo(i,2,N){
if (!bz[i]){
pri[++top]=i;
mu[i]=-1;
}
fo(j,1,top){
if ((ll)i*pri[j]>N) break;
bz[i*pri[j]]=1;
if (i%pri[j]==0){
mu[i*pri[j]]=0;
break;
}
mu[i*pri[j]]=-mu[i];
}
}
fo(i,1,N)
fo(j,1,N/i)
d[i*j]++;
fo(i,1,N) sum[i]=sum[i-1]+d[i],mem[i]=-1;
}
int calcg(ll x){
if (x<=N) return sum[x];
if (mem[n/x]>=0) return mem[n/x];
ll i=1,j;
int t=0;
while (i<=x){
j=x/(x/i);
(t+=(ll)(j-i+1)*(x/i)%mo)%=mo;
i=j+1;
}
return mem[n/x]=t;
}
int calcf(ll x){
ll i=1,j;
int t=0;
while (i<=x){
j=x/(x/i);
(t+=(ll)calcg(x/i)*(j-i+1)%mo)%=mo;
i=j+1;
}
return t;
}
int main(){
freopen("ra.in","r",stdin);freopen("ra.out","w",stdout);
prepare();
scanf("%lld",&n);
fo(i,1,floor(sqrt(n)))
if (mu[i]) (ans+=mu[i]*calcf(n/i/i))%=mo;
ans=((ll)(n%mo)*(n%mo)%mo-ans)%mo;
(ans+=mo)%=mo;
printf("%d\n",ans);
}