[51NOD1237]最大公约数之和 V3

题目大意

给定$n$，试求

$$\sum_{i=1}^n\sum_{j=1}^n\gcd(i,j)$$
结果对$10^9+7$取模。

$2\le n\le10^{10}$

---

题目分析

我们将题目改为求

$$\sum_{i=1}^n\sum_{j=1}^i\gcd(i,j)$$
然后将答案乘二再减去$1$到$n$的和即可。
那么上面这条式子是什么呢？
$$\begin{align} \sum_{i=1}^n\sum_{j=1}^i\gcd(i,j)&=\sum_{d=1}^nd\sum_{i=1}^{\lfloor\frac nd\rfloor}\varphi(i)\\ &=\sum_{d=1}^n\varphi(d){\left(1+\lfloor\frac nd\rfloor\right)\lfloor\frac nd\rfloor\over2} \end{align}$$
对$\lfloor\frac nd\rfloor$分块，对前面的欧拉函数使用杜教筛求前缀和。
时间复杂度？这个算起来有点麻烦，我来口胡一下。
假设我们预处理的部分大小为$a$
$$\sum_{i=1}^\sqrt n\sum_{j=1}^{\lfloor\frac ia\rfloor}\mathrm O\left(\sqrt{\frac ij}\right)+\sum_{i=1}^\sqrt n\sum_{j=1}^{\lfloor\frac n{ia}\rfloor}\mathrm O\left(\sqrt\frac n{ij}\right)$$
显然右边部分的复杂度要大于左边，使用积分来计算：
$$\begin{align} \sum_{i=1}^\sqrt n\sum_{j=1}^{\lfloor\frac n{ia}\rfloor}\mathrm O\left(\sqrt\frac n{ij}\right)&\approx\mathrm O\left(\int_0^\sqrt n\left(\int_0^\frac n{ax}\sqrt\frac n{xy}\ \mathrm dy\right)\ \mathrm dx\right)\\ &=\mathrm O\left(\int_0^\sqrt n\frac nx\sqrt\frac1a\right)\\ &=\mathrm O\left(\sqrt\frac1an\log n\right) \end{align}$$
然后$a$可以在空间限制内尽量取大一点，我取了$1.5\times10^7$。

---

代码实现

我的代码莫名其妙常数巨大，卡着空间取的让我卡着时限过了这题。。。

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

using namespace std;

typedef long long LL;

const int P=1000000007;
const int itwo=500000004;
const int L=15000000;
const int N=100050;

int pri[L+5],phi[L+5],f[L+5],sum[L+5];
int S[N],vis[N];
bool mark[L+5];
int ans,l,s,cnt;
LL n;

void pre()
{
    s=trunc(sqrt(n)),l=L;
    mark[1]=1,sum[1]=phi[1]=1,f[1]=1;
    for (int i=2;i<=l;++i)
    {
        if (!mark[i]) pri[++pri[0]]=i,phi[i]=i-1,f[i]=i;
        for (int j=1,k;j<=pri[0];++j)
        {
            if (1ll*pri[j]*i>l) break;
            mark[k=pri[j]*i]=1,f[k]=pri[j];
            phi[k]=phi[i]*(pri[j]-(f[i]!=f[k]));
            if (!(i%pri[j])) break;
        }
        sum[i]=(sum[i-1]+phi[i])%P;
    }
}

int id(LL x){return n/x;}

int sieve(LL n)
{
    if (n<=l) return sum[n];
    int idn=id(n);
    if (vis[idn]==cnt) return S[idn];
    int res=1ll*((n+1)%P)*(n%P)%P*itwo%P;
    for (LL st=2,en,x;st<n;st=en)
    {
        x=n/st,en=n/x+1;
        res=(res-1ll*sieve(x)*(en-st)%P+P)%P;
    }
    return vis[idn]=cnt,S[idn]=res;
}

int main()
{
    freopen("gcdsumIII.in","r",stdin),freopen("gcdsumIII.out","w",stdout);
    scanf("%lld",&n),pre();
    ans=0;
    for (LL st=1,en,x;st<n;st=en)
    {
        x=n/st,en=n/x+1;
        ans=(ans+1ll*((x+1)%P)%P*(x%P)%P*itwo%P*(((++cnt,sieve(en-1))-(++cnt,sieve(st-1))+P)%P)%P)%P;
    }
    ans=(ans*2%P-1ll*((n+1)%P)*(n%P)%P*itwo%P+P)%P;
    printf("%d\n",ans);
    fclose(stdin),fclose(stdout);
    return 0;
}