[JZOJ4872]太阳神

题目大意

求满足如下条件的数对$(a,b)$对数：$a,b\in \mathbb Z^+$且$a,b\le n$而$\mathrm{lcm}(a,b)>n$。
答案对$10^9+7$取模。

$1\le n\le 10^{10}$

---

题目分析

正难则反：用$n^2$减去$\mathrm{lcm}(a,b)\le n$的对数。
我们枚举$d=\gcd(a,b)$，令$a'={a\over d},b'={b\over d}$。
显然$a'b'd\le n$。即我们要求的是$\gcd(a',b')=1$且$a'b'd\le n$的三元组$(a',b',d)$个数。

$$ans=\sum_{d=1}^n\sum_{a'=1}^{\lfloor{n\over d}\rfloor}\sum_{b'=1}^{\lfloor{n\over da'}\rfloor}[\gcd(a',b')=1]$$
考虑莫比乌斯反演，枚举 $i=\gcd(a',b')$，令 $a''={a'\over i},b''={b'\over i}$：
$$ans=\sum_{d=1}^n\sum_i\mu(i)\sum_{a''=1}\sum_{b''=1}\lfloor{n\over da''i^2}\rfloor$$
这个看起来很难看，我们更换主体：
$$ans=\sum_{i=1}^\sqrt n\mu(i)[a''b''d\le{n\over i^2}的三元组(a'',b'',d)个数]$$
这里 $i$是小于等于$\sqrt n$的，因此我们可以线筛出 $\mu$函数的值。
现在问题在于如何求满足 $abc\le m$的三元组个数。
我们可以限定 $a\le b\le c$，先计算满足这个条件的三元组个数。显然 $a$可以在$m^{1\over3}$的范围内枚举， $b$可以在$\sqrt{m\over a}$的范围内枚举， $c$的个数可以直接$\mathrm O(1)$计算。这样计算一次的复杂度是多少呢？使用定积分：
$$\begin{align} &\int_1^{m^{1\over3}}\sqrt{m\over x}dx\\ =&\sqrt m\int_1^{m^{1\over3}}x^{-{1\over2}}dx\\ =&\sqrt m\left(2(m^{1\over3})^{1\over2}-2\right)\\ =&\sqrt m\mathrm O(m^{1\over6})\\ =&\mathrm O(m^{2\over3}) \end{align}$$
但是注意到我们还要枚举 $i$，这样时间复杂度会不会特别大呢？还是使用定积分：
$$\begin{align} &\int_1^\sqrt n\left({n\over x^2}\right)^{2\over3}dx\\ =&n^{2\over3}\int_1^\sqrt nx^{-{4\over3}}dx\\ =&n^{2\over3}\left(-3n^{-{1\over6}}+3\right)\\ =&\mathrm O(n^{2\over3})（上式后面部分是一个小于3的常数） \end{align}$$
但是这样并不是所有的方案数。因此我们要将其乘上 $3!$，但是你发现这样你会算重。于是你可以枚举两个相等的情况（ $\mathrm O(n^{1\over2})$的时间复杂度显然是小于 $\mathrm O(n^{2\over3})$的），然后减一减，再减去三个相等的情况。具体细节就交给你们自己讨论了。
总的时间复杂度 $\mathrm O(n^{2\over3})$。

---

代码实现

#include <iostream>
#include <cstdio>
#include <cmath>

using namespace std;

typedef long long LL;

const int P=1000000007;
const int A=100000;

int mu[A+50],pri[A+50],f[A+50];
int ans;
LL n;

void calc()
{
    mu[1]=1;
    for (int i=2;i<=A;i++)
    {
        if (!f[i]) mu[i]=-1,f[pri[++pri[0]]=i]=i;
        for (int j=1;j<=pri[0];j++)
        {
            if (1ll*i*pri[j]>A) break;
            f[i*pri[j]]=pri[j];
            if (f[i]==pri[j]) mu[i*pri[j]]=0;
            else mu[i*pri[j]]=-mu[i];
            if (!(i%pri[j])) break;
        }
    }
    int limi=trunc(sqrt(n));
    for (int i=1,cnt,lima;i<=limi;i++)
    {
        if (!mu[i]) continue;
        LL x=n/(1ll*i*i);
        cnt=0;
        for (int a=1,limb;1ll*a*a*a<=x;a++)
        {
            limb=trunc(sqrt(x/a));
            (((cnt-=2)%=P)+=P)%=P;
            for (int b=a;b<=limb;b++) (cnt+=6ll*(x/(1ll*a*b)-b+1)%P)%=P;
        }
        lima=trunc(sqrt(x));
        for (int a=1;a<=lima;a++) (((cnt-=1ll*(x/(1ll*a*a))*3%P)%=P)+=P)%=P;
        (((ans+=1ll*cnt*mu[i])%=P)+=P)%=P;
    }
    ans=((1ll*(n%P)*(n%P)%P-ans)%P+P)%P;
}

int main()
{
    freopen("ra.in","r",stdin),freopen("ra.out","w",stdout);
    scanf("%lld",&n);
    calc();
    printf("%d\n",ans);
    fclose(stdin),fclose(stdout);
    return 0;
}