51Nod1237：最大公约数之和

51Nod1237：最大公约数之和

题意：

• 求\(\sum_{i=1}^N\sum_{j=1}^Ngcd(i,j)\)。
• 数据范围\(:N\leq 10^{10}\)。

思路：

• 首先先来推一下\(\sum_{i=1}^N\sum_{j=1}^N gcd(i,j)\)的公式。

• 原式=\(\sum_{d=1}^Nd\sum_{i=1}^N\sum_{j=1}^N[gcd(i,j)==d]\)。

• \(=\sum_{d=1}^Nd\sum_{i=1}^{\frac{N}{d}}\sum_{j=1}^{\frac{N}{d}}[gcd(i,j)=1]\)。

• \(=\sum_{d=1}^nd\sum_{i=1}^{\frac{N}{d}}2\phi(i)-1\)。

• 设\(f(i)=\sum_{1}^i\phi(i)\)

• 则原式有：

• \(\sum_{i=1}^ni(2f(\frac{n}{i})-1\))。

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int maxn = 5e6;
const ll mod = 1e9 + 7;
const int inv2 = (mod+1)/2;
bool vis[maxn+10];
int primes[maxn+10], cnt;
ll phi[maxn+10];
ll n;
void init(int n)
{
    phi[1] = 1;
    for(int i = 2; i <= n; i++)
    {
        if(!vis[i])
        {
            primes[++cnt] = i;
            phi[i] = (i-1);
        }
        for(int j = 1; primes[j] <= n/i; j++)
        {
            vis[primes[j]*i] = 1;
            if(i % primes[j] == 0)
            {
                phi[primes[j]*i] = (phi[i]*primes[j])%mod;
                break;
            }
            else phi[primes[j]*i] = (phi[i]*(primes[j]-1))%mod;
        }
    }
    for(int i = 1; i <= n; i++)
        phi[i] += phi[i-1], phi[i] %= mod;
}

unordered_map<ll, ll> Sphi;

inline int getSphi(ll n)
{
    if(n <= maxn) return phi[n];
    if(Sphi[n]) return Sphi[n];
    ll res = 1ll*n%mod*((n+1)%mod)%mod*inv2%mod;
    for(ll l = 2,r; l <= n; l = r+1)
    {
        r = n/(n/l);
        res -= 1ll*(r-l+1)%mod*getSphi(n/l)%mod;
        res = (res+mod)%mod;
    }
    return Sphi[n] = res;
}

int main()
{
    scanf("%lld", &n);
    init(maxn);
    ll ans = 0;
    for(ll l = 1, r; l <= n; l = r + 1)
    {
        r = n/(n/l);
        ans=(ans+(1ll*(r-l+1)%mod*((r+l)%mod)%mod*inv2%mod)%mod*((2ll*getSphi(n/l)%mod)-1+mod)%mod)%mod;
    }
    printf("%lld\n", ans);
    return 0;
}