51nod1237(EES解法，省空间)

这里用的是最近比较流行的EES方法解决的
不了解的可以看下这篇文章：
https://blog.csdn.net/Feynman1999/article/details/82874491
这种方法只要能找出$f(p^e)$的表达式，并且当$e+1$时可以$O(1)$维护，就可以用EES筛法去做
对于本题而言，考虑函数$S(n)={\sum }_{i=1}^{n}gcd(i,n)$，可知若能快速求出其和函数的值记做res，则$2\ast res-(1+2+3+...+n)$即为ans
如何求res呢？
可知函数$S(n)$是积性函数，实际上$S(n)={\sum }_{d|N}\phi (\frac{N}{d})\ast d$，可以看到他是欧拉函数和单位函数的狄利克雷卷积。那我们列出$S(p)=2p-1$,$S(p^{e+1})=p\ast S(P^e)+\phi (p^{e+1})$ ，那只要维护一个欧拉函数即可，$\phi (p)=p-1,\phi (p^{e+1})=p\ast \phi (p^e)$
维护$p^1,p^0$的前缀即可（看了EES才知道这里前缀是什么意思）
没有特别要注意的地方

代码

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int mod=1e9+7;
const int ni6=166666668;
const int ni2=500000004;
ll n,M;
vector<int> pre[2],hou[2],primes,pri_ni;

inline int add(const int x, const int v) {
    return x + v >= mod ? x + v - mod : x + v;
}
inline int dec(const int x, const int v) {
    return x - v < 0 ? x - v + mod : x - v;
}

//这里res是n/枚举的数
int dfs(ll res, int last, ll f){
    //最大质因子是prime[last-1] 但将1放在外面值显然一样
    //cout<<n<<" "<<res<<endl;
    int t=dec((res > M ? hou[1][n/res] : pre[1][res]),pre[1][primes[last]-1]);
    //cout<<t<<endl;
    int ret= (ll)t*f%mod;//这里需修改
    for(int i=last;i<(int) primes.size();++i){
        int p = primes[i];
        if((ll)p*p > res) break;
        ll phi=p-1;
        ll ftp=(2*p-1);
        for(ll q=p,nres=res,nf=f*ftp%mod;q*p<=res;q*=p){//nf需修改
            ret = add(ret,dfs(nres/=p,i+1,nf));//枚举更大的数
            phi = phi*p%mod;
            ftp = add(p*ftp%mod,phi);
            nf = f*ftp%mod;//继续枚举当前素数的指数
            ret =add(ret,nf);//指数大于1时，记上贡献
        }
    }
    return ret;
}

inline int ff(ll x){
    x%=mod;
    return x*(x+1)%mod*ni2%mod;
}

inline int fff(ll x){
    x%=mod;
    return x*(x+1)%mod*(2*x+1)%mod*ni6%mod;
}

int solve(ll n){
    M=sqrt(n);
    for(int i=0;i<2;++i){
        pre[i].clear();pre[i].resize(M+1);
        hou[i].clear();hou[i].resize(M+1);
    }
    primes.clear();primes.reserve(M+1);
    pri_ni.clear();pri_ni.reserve(M+1);
    for(int i=1;i<=M;++i){
        pre[0][i]=i-1;
        hou[0][i]=(n/i-1)%mod;
        pre[1][i]=dec(ff(i),1);;
        hou[1][i]=dec(ff(n/i),1);
    }
    for(int p=2;p<=M;++p){
        if(pre[0][p]==pre[0][p-1]) continue;
        primes.push_back(p);
        const ll q=(ll)p*p,m=n/p;
        const int pnt0=pre[0][p-1],pnt1=pre[1][p-1];
        const int mid=M/p;
        const int End=min((ll)M,n/q);
        for(int i=1;i<=mid;++i){
            hou[0][i]=dec(hou[0][i],dec(hou[0][i*p],pnt0));
            hou[1][i]=dec(hou[1][i],dec(hou[1][i*p],pnt1)*(ll)p%mod);
            //hou[2][i]=dec(hou[2][i],dec(hou[2][i*p],pnt2)*q%mod);
        }
        for(int i=mid+1;i<=End;++i){
            hou[0][i]=dec(hou[0][i],dec(pre[0][m/i],pnt0));
            hou[1][i]=dec(hou[1][i],dec(pre[1][m/i],pnt1)*(ll)p%mod);
            //hou[2][i]=dec(hou[2][i],dec(pre[2][m/i],pnt2)*q%mod);
        }
        for(int i=M;i>=q;--i){
            pre[0][i]=dec(pre[0][i],dec(pre[0][i/p],pnt0));
            pre[1][i]=dec(pre[1][i],dec(pre[1][i/p],pnt1)*(ll)p%mod);
            //pre[2][i]=dec(pre[2][i],dec(pre[2][i/p],pnt2)*q%mod);
        }
    }
    primes.push_back(M+1);
    for (int i = 1; i <= M; i++) {
        pre[1][i] = dec(2*pre[1][i]%mod, pre[0][i]);//2*p-1
        hou[1][i] = dec(2*hou[1][i]%mod, hou[0][i]);
    }
    return n>1 ? add(dfs(n,0,1),1) : 1;
}

int main()
{
    //freopen("in.txt","r",stdin);
    ios::sync_with_stdio(false);
    cin>>n;
    cout<<((2LL*solve(n)-ff(n))%mod+mod)%mod<<endl;
    return 0;
}

另外附上杜教筛的代码：

#include<bits/stdc++.h>
#define pb push_back
#define mp make_pair
#define bit(a,b) ((a>>b)&1) //from 0
#define all(x) (x).begin(),(x).end()
const int ni2=500000004;
using namespace std;
typedef long long ll;
int debug_num=0;

const int mod=1e9+7;
const int maxn=1e7+10;
bool valid[maxn];
int ans[maxn/10];
ll phi[maxn];
const int maxm=(1LL<<37)/maxn;
ll help1[maxm];
bool vis[maxm];
ll tot,up,m;

inline ll add(const ll x, const ll v) {
    //return x + v >= mod ? x + v - mod : x + v;
    return (x+v)%mod;
}

inline ll dec(const ll x, const ll v) {
    //return x - v < 0 ? x - v + mod : x - v;
    return (x-v)%mod;
}

void get_prime(int n)
{
    memset(valid,true,sizeof(valid));
    tot=0;
    phi[1]=1;
    for(int i=2;i<=n;++i){
        if(valid[i]){
            ans[++tot]=i;
            phi[i]=i-1;
        }
        for(int j=1;j<=tot && ans[j]*i<=n;++j){
            int tp=ans[j]*i;
            valid[tp]=false;
            if(i%ans[j]==0){
                phi[tp]=phi[i]*ans[j];
                break;
            }
            else{
                phi[tp]=phi[i]*(ans[j]-1);
            }
        }
    }
    for(int i=1;i<=n;++i){
        phi[i]=add(phi[i],phi[i-1]);
    }
}

ll get_phi(ll n)
{
    return (n<=up) ? phi[n] : help1[m/n] ;
}

inline  ff(ll x){
    x%=mod;
    return x*(x+1)%mod*ni2%mod;
}

void solve1(ll n)
{
    int t;
    if(n<=up || vis[t=m/n]) return ;
    vis[t]=true;
    help1[t]=ff(n);
    for(ll l=2,r;l<=n;l=r+1){
        r=n/(n/l);
        solve1(n/r);
        help1[t]=dec(help1[t]%mod,(r-l+1)%mod*get_phi(n/r)%mod);
    }
}


ll solve2(ll n)
{
    ll ans=0;
    for(ll l=1,r;l<=n;l=r+1){
        r=n/(n/l);
        ll tp=n/r;
        if(tp<=up){
            ans=add(ans,(ff(r)-ff(l-1))*phi[tp]%mod);
        }
        else{
            ans=add(ans,(ff(r)-ff(l-1))*help1[m/tp]%mod);
        }
    }
    return ans;
}

int main()
{
    //freopen("in.txt","r",stdin);
    up=maxn-10;
    get_prime(up);
    //cout<<clock()<<endl;
    ll n;
    cin>>n;
    m=n;
    if(m>up){
        memset(vis,false,sizeof(vis));
        solve1(n);
    }
    cout<<((2*solve2(n)%mod-ff(n))%mod+mod)%mod<<endl;
    return 0;
}