【51nod 1237】 最大公约数之和 V3

description

给出一个数N（<= 10^10），输出小于等于N的所有数，两两之间的最大公约数之和。

相当于计算这段程序（程序中的gcd(i,j)表示i与j的最大公约数）：
由于结果很大，输出Mod 1000000007的结果。

G=0;
for(i=1;i< N;i++)
for(j=1;j<=N;j++)
{
G = (G + gcd(i,j)) % 1000000007;
}

Solution

这道题想出来后被人告知是杜教筛……
前面的都很顺利，非常快的推出式子。

$$ans=\sum_{i=1}^{n}\sum_{j=1}^{n}gcd(i,j)=\sum_{d=1}^{n}d\sum_{i=1}^{n/d}{{\lfloor{n\over id }\rfloor}*{\lfloor{n\over id }\rfloor}*\mu(i)}=\sum_{T=1}^{n}{\lfloor{n\over T }\rfloor}*{\lfloor{n\over T }\rfloor}*\sum_{d|T}d*\mu(T/d)$$
推到这里，我们发现一个很喜欢的式子 $\sum_{d|T}d*\mu(T/d)$，因为我们知道 $\phi(n)=\sum_{d|n}d*\mu(n/d)$。
$$ans=\sum_{T=1}^{n}{\lfloor{n\over T }\rfloor}*{\lfloor{n\over T }\rfloor}*\phi(T)。$$
然后考虑分块，所以我们要求一段数的欧拉函数，这和之前做的一道 51nod1239 欧拉函数之和简直一模一样，这里就不再推导了，直接上结论。
$设f(x)=\sum_{i=1}^{x}\phi(i)$
那么 $$f(x)=(1+x)*x/2-\sum_{T=2}^{x}f(\lfloor x/T \rfloor)。$$ 所以ans中欧拉函数再分块中用 $f(j)-f(i-1)$算一下即可。复杂度O( $n^{2/3}$)。

Code

#include<iostream>
#include<cmath>
#include<cstring>
#include<cstdio>
#include<algorithm>
#define ll long long
using namespace std;
const ll maxn=1e6+5,mo=1e9+7,mo2=5e8+4;
ll d[maxn],bz[maxn],p[maxn],h[maxn*10],f[maxn*10];
ll n,i,t,j,k,l,x,y,z,ans;
int hash(ll x){
    ll t=x%maxn;
    while (h[t] && h[t]!=x) t=(t+1)%maxn;
    return t;
}
ll dg(ll n){
    if (n<=maxn) return p[n];
    ll t=n%mo,k=t*(t+1)%mo*mo2%mo,i=2,x=hash(n),y;
    if (h[x]) return f[x];
    while (i<=n){
        t=n/(n/i);
        k-=dg(n/i)*(t-i+1)%mo;i=t+1;
    }
    k=(k%mo+mo)%mo;h[x]=n;f[x]=k;return k;
}
int main(){
    //freopen("data.in","r",stdin);
    scanf("%lld",&n);p[1]=1;
    for (i=2;i<=maxn;i++){
        if (!bz[i]) d[++d[0]]=i,p[i]=i-1;
        for (j=1;j<=d[0];j++){
            if (i*d[j]>maxn) break;
            bz[i*d[j]]=1;
            if (i%d[j]==0){
                p[i*d[j]]=p[i]*d[j];break;
            }else p[i*d[j]]=p[i]*p[d[j]];
        }
    }
    for (i=1;i<=maxn;i++)
        p[i]+=p[i-1];
    i=1;
    while (i<=n){
        t=n/(n/i);
        k=(dg(t)-dg(i-1)+mo)%mo;l=n/i%mo;
        ans=(ans+l*l%mo*k%mo)%mo;
        i=t+1;
    }
    printf("%lld\n",ans);
}