[莫比乌斯反演+杜教筛] 51Nod1237: 最大公约数之和 V3

题意

求$\sum_i^n\sum_j^ngcd(i,j)$
$n\le10^{10}$

题解

$$\sum_i^n\sum_j^ngcd(i,j)=\sum_{d=1}^n\sum_{i=1}^n\sum_{j=1}^n[gcd(i,j)=d]*d=\sum_{d=1}^n\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{d}\rfloor} [gcd(i,j)=1]*d$$
然后反演：
$$=\sum_{d=1}^n\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{n}{d}\rfloor} \sum_{k|gcd(i,j)}\mu(k)*d=\sum_{k=1}^n\mu(k)\sum_{d=1}^n\sum_{k|i}^{\lfloor\frac{n}{d}\rfloor}\sum_{k|j}^{\lfloor\frac{n}{d}\rfloor}d=\sum_{k=1}^n\mu(k)\sum_{d=1}^n \lfloor\frac{n}{dk}\rfloor \lfloor\frac{n}{dk}\rfloor d$$
设 $T=dk$,得
$$=\sum_{T=1}^n \lfloor\frac{n}{T}\rfloor \lfloor\frac{n}{T}\rfloor \sum_{d|T} d*\mu(\frac{T}{d})=\sum_{T=1}^n \lfloor\frac{n}{T}\rfloor \lfloor\frac{n}{T}\rfloor \phi(T)$$
然后就好求了，分块求和即可。
还要用杜教筛求 $\phi(i)$的前缀和，具体见 此题。

有个地方没模到wa了无数次……这种问题要小心，大数据才看的出来。
不过自己手推感觉还是爽的。

#include<cstdio>
#include<cstring>
#include<algorithm>
#include<tr1/unordered_map>
using namespace std;
using namespace std::tr1;
typedef long long LL;
const LL N=10000000, maxn=N+5, MOD=1000000007, inv=500000004;
int ans,p[maxn],phi[maxn];
LL n; 
bool vis[maxn];
void get_phi(){
    memset(vis,1,sizeof(vis)); 
    phi[1]=1;
    for(int i=2;i<=N;i++){
        if(vis[i]) p[++p[0]]=i, phi[i]=i-1;
        for(int j=1;j<=p[0]&&(LL)i*p[j]<=N;j++){
            vis[i*p[j]]=false;
            if(i%p[j]==0){ phi[i*p[j]]=phi[i]*p[j]; break; }
            phi[i*p[j]]=phi[i]*phi[p[j]];
        }
    }
    for(int i=1;i<=N;i++) (phi[i]+=phi[i-1])%=MOD;
}
unordered_map< LL, int > lst;
int Sum_phi(LL n){
    if(n<=N) return phi[n];
    if(lst.find(n)!=lst.end()) return lst[n];
    int res=(LL)n%MOD*((n+1)%MOD)%MOD*inv%MOD;
    for(LL i=2,nxt;i<=n;i=nxt+1){
        nxt=n/(n/i);
        (res+=MOD-(LL)(nxt-i+1)*Sum_phi(n/i)%MOD)%=MOD;
    }
    return lst[n]=res;
}
int main(){
    freopen("51nod1237.in","r",stdin);
    freopen("51nod1237.out","w",stdout);
    get_phi();
    scanf("%lld",&n);
    for(LL i=1,nxt;i<=n;i=nxt+1){
        nxt=n/(n/i); 
        (ans+=((n/i)%MOD)*((n/i)%MOD)%MOD*(Sum_phi(nxt)+MOD-Sum_phi(i-1))%MOD)%=MOD; 
    }
    printf("%d\n",ans);
    return 0;
}