【51nod1237】 最大公约数之和 V3 (杜教筛)

基准时间限制：5 秒 空间限制：262144 KB 分值: 640 难度：8级算法题
给出一个数N，输出小于等于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;
}
Input
输入一个数N。(2 <= N <= 10^10)
Output
输出G Mod 1000000007的结果。
Input示例
100
Output示例
31080

这个其实要比最小公倍数的版本要好做，之前我也推过，很容易我们可以得到：

$$ans=\sum_{d=1}^n\lfloor\frac{n}{d}\rfloor^2\phi(d)$$
现在只是n太大了，那么在这里我们仍然可以分块来处理，只是需要处理欧拉函数的前缀和，那这个用杜教筛很容易就解决了。有： $$sum(n)=\frac{n^2+n}{2}-\sum_{i=2}^nsum(\lfloor\frac{n}{k}\rfloor)$$
代码：

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<map>
#include<cmath>
#define maxx 5000000
#define mod 1000000007
#define ll long long
using namespace std;
bool isP[maxx];
int prime[400000];
int cnt;
ll phi[maxx];
ll inv2=500000004;
void init()
{
    phi[1]=1;
    for(int i=2;i<maxx;i++)
    {
        if(!isP[i])
        {
            prime[cnt++]=i;
            phi[i]=i-1;
        }
        for(int j=0;j<cnt&&(ll)i*prime[j]<maxx;j++)
        {
            isP[i*prime[j]]=true;
            if(i%prime[j])
                phi[i*prime[j]]=phi[i]*(prime[j]-1);
            else
            {
                phi[i*prime[j]]=phi[i]*prime[j];
                break;
            }
        }
    }
    for(int i=1;i<maxx;i++)
        phi[i]=(phi[i]+phi[i-1])%mod;
}
map<ll,ll>M;
ll work(ll x)
{
    if(x<maxx)return phi[x];
    if(M[x])return M[x];
    ll t=x%mod;
    ll ans=t*(t+1)%mod*inv2%mod;
    for(ll i=2,last;i<=x;i=last+1)
    {
        last=x/(x/i);
        ans=(ans-(last-i+1)%mod*work(x/i)%mod)%mod;
    }
    if(ans<0)ans=(ans+mod)%mod;
    M[x]=ans;
    return ans;
}
int main()
{
    init();
    ll ans=0;
    ll n;
    cin>>n;
    for(ll i=1,last;i<=n;i=last+1)
    {
        last=n/(n/i);
        ll now=(n/i)%mod;
        now=now*now%mod;
        ll base=(work(last)-work(i-1)+mod)%mod;
        ans=(ans+now*base%mod)%mod;
    }
    cout<<ans<<endl;
    return 0;
}