[日常摸鱼]51nod1237-最大公约数之和V3-杜教筛

题意：求$\sum_{i=1}^n \sum_{j=1}^n gcd(i,j),n<=1e10$

 

---

 

 

之前刚好在UVA上也做过一个这样求和的题目，不过那个数据范围比较小，一开始用类似的方法

$ans=\sum_{i=1}^n \sum_{j=1}^i gcd(i,j)-\sum_{i=1}^n i$

先考虑化简$\sum_{i=1}^n gcd(i,n)$变成好求和的形式

$$\begin{aligned} \sum_{i=1}^n gcd(i,n) &=\sum_{i=1}^n \sum_{d=1}^n d*[gcd(i,n)=d]\\ &=\sum_{d=1}^n d \sum_{i=1}^n [gcd(i,n)=d] \\ &=\sum_{d=1}^n d \sum_{\frac{i}{d}=1}^{\frac{n}{d}} [gcd(\frac{i}{d},\frac{n}{d})=1]\\ &=\sum_{d|n} d*\phi(\frac{n}{d}) \end{aligned}$$

发现是$f(n)=n$和$g(n)=\phi(n)$的卷积，令$S(n)$表示$g$的前缀和，然后非常套路地，然后刚刚那一坨的前缀和就变成求$\sum_{i=1}^n i*S(\lfloor \frac{n}{i} \rfloor)$，欧拉函数前缀和直接用杜教筛算，然后分块求和

 

一开始取模一直写挂…orz

（有点懒直接用map存了）

#include<cstdio>
#include<cstring>
#include<map>
using namespace std;

typedef long long lint;

const lint MOD=1000000007;
const lint N=4000005;
const lint G=100005;

lint n,tot,inv2;
lint pri[N],phi[N];
bool p[N];
map<lint,lint>mp;

inline void init()
{
    p[1]=1;phi[1]=1;
    for(register lint i=2;i<N;i++)
    {
        if(!p[i])
        {
            pri[++tot]=i;
            phi[i]=i-1;
        }
        for(register lint j=1;j<=tot&&i*pri[j]<N;j++)
        {
            lint t=i*pri[j];p[t]=1;
            if(i%pri[j]==0){phi[t]=phi[i]*pri[j];break;}
            phi[t]=phi[i]*(pri[j]-1);
        }
    }
    for(register lint i=1;i<N;i++)phi[i]=(phi[i]+phi[i-1])%MOD;
}

inline lint pow_mod(lint a,lint b,lint p)
{
    lint res=1;
    for(;b;b>>=1,a=(a*a)%p)if(b&1)res=(res*a)%p;
    return res%p;
}

inline lint sum(lint x)
{
    return x%MOD*((x+1)%MOD)%MOD*inv2%MOD;
}

inline lint calc_phi(lint x)
{
    if(x<N)return phi[x];
    if(mp.count(x))return mp[x];
    lint res=sum(x),pos;
    for(register lint i=2;i<=x;i=pos+1)
    {
        pos=x/(x/i);
        res-=((pos-i+1)%MOD*calc_phi(x/i)%MOD)%MOD;
        res=(res%MOD+MOD)%MOD;
    }return mp[x]=res;
}

inline lint calc_ans(lint x)
{
    lint res=0,pos;
    for(register lint i=1;i<=x;i=pos+1)
    {
        pos=x/(x/i);
        res+=(sum(pos)-sum(i-1))%MOD*calc_phi(x/i)%MOD;
        res=(res%MOD+MOD)%MOD;
    }return res;
}

int main()
{
    init();inv2=pow_mod(2,MOD-2,MOD);
    scanf("%lld",&n);lint ans=calc_ans(n)%MOD;
    ans=(ans*2)%MOD-sum(n);ans=(ans%MOD+MOD)%MOD;
    printf("%lld",ans);
    return 0;
}