太阳神

题目描述

太阳神拉很喜欢最小公倍数，有一天他想到了一个关于最小公倍数的题目。
求满足如下条件的数对(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的取值就有$\frac{n}{ab}$
$ans=\sum_{a=1}^n\sum_{b=1}^n[(a,b)=1]*\frac{n}{ab}$
我们尝试莫比乌斯反演一下，也就是
设$f[k]=\sum_{a=1}^n\sum_{b=1}^n[(a,b)=k]*\frac{n}{ab}$
$g[k]=\sum_{i=1}^{\frac{n}{k}}f[k*i]$
反演一下变成
$f[k]=\sum_{i=1}^n\mu(i)*g[k*i]$
g是很容易算的，所以
$f[k]=\sum_{i=1}^n\mu(i)*\sum_{a'=1}^{\frac{n}{k*i}}\sum_{b'=1}^{\frac{n}{k*i}}\frac{(\frac{n}{k^2})}{a'b'}$
现在我们需要f[1]，代入k=1得
$ans=\sum_{d=1}^n\mu(d)*\sum_{a'=1}^{\frac{n}{d^2}}\sum_{b'=1}^{\frac{n}{d^2}}\frac{(\frac{n}{d^2})}{a'b'}$
观察后面，可以改变d的枚举上界
$ans=\sum_{d=1}^{\sqrt n}\mu(d)*\sum_{a'=1}^{\frac{n}{d^2}}\sum_{b'=1}^{\frac{n}{d^2}}\frac{(\frac{n}{d^2})}{a'b'}$
后面那个东西只与$\frac{n}{d^2}$有关，设其为$f(\frac{n}{d^2})$

算f

设$f(x)=\sum_{i=1}^x\sum_{j=1}^x\frac{x}{ij}$
可以看出可以根据$\frac{x}{i}$的取值分块！每一块答案显然相同，块数为$\sqrt x$
假如i确定了，那么答案为$\sum_{j=1}^{\frac{x}{i}}\frac{(\frac{x}{i})}{j}$
后面那个东西只与$\frac{x}{i}$有关，设其为$g(\frac{x}{i})$

算g

设$g(x)=\sum_{i=1}^x\frac{x}{i}$
我们应该怎么算g呢？
有这样一个奇妙的方法。
首先，先知道一个结论。
$g(x)=\sum_{i=1}^xd(i)$
d函数表示约数个数。
为什么呢？你考虑$\frac{x}{i}$，那么说明$1*i,2*i,3*i……\frac{x}{i}*i$都在x内，而且i是上述提到的$\frac{x}{i}$个数的约数。
因此可以表示为约数个数前缀和。
那么也就是说g可以线筛？
但显然不能全都线筛啊？
dalao说，以$n^{2/3}$为界，记为N。我们需要筛出g(1~N)，然后对于x<=N便直接得解。
如果x>N？观察g的算式，我们显然可以按照取值分块，得到一个根号算法。那我们就这样做最快咯。不过这样会T，加个记忆化即可。
怎么记忆化呢？dalao说，x可以映射到n/x，这样一定不会重复。这样就可以开小记忆化数组。
总复杂度？dalao说，是$O(n^{2/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);
}