51NOD 1237 最大公约数之和 V3 [杜教筛]

1237 最大公约数之和 V3

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

---

令\(A(n)=\sum_{i=1}^n(n,i) = \sum_{d\mid n}d \cdot \varphi(\frac{n}{d})\)

\(ans = 2*\sum_{i=1}^n A(i) -\sum_{i=1}^ni\)

套路推♂倒
\[ S(n) =\sum_{i=1}^n\sum_{d\mid i}d \cdot \varphi(\frac{i}{d}) =\sum_{i=1}^n i \sum_{d=1}^{\lfloor \frac{n}{i} \rfloor} \varphi(d) \]
杜教筛\(\varphi\)的前缀和后整除分块

#include <iostream>
#include <cstdio>
#include <cstring>
#include <algorithm>
#include <cmath>
using namespace std;
typedef long long ll;
const int N = 4641590, U = 4641588, mo = 1e9+7, inv2 = 500000004;
inline ll read(){
    char c=getchar(); ll x=0,f=1;
    while(c<'0' || c>'9') {if(c=='-')f=-1; c=getchar();}
    while(c>='0' && c<='9') {x=x*10+c-'0'; c=getchar();}
    return x*f;
}

ll n;
inline void mod(ll &x) {if(x>=mo) x-=mo; else if(x<0) x+=mo;}
bool notp[N]; int p[N/10]; ll phi[N];
void sieve(int n) {
    phi[1]=1; 
    for(int i=2; i<=n; i++) {
        if(!notp[i]) p[++p[0]] = i, phi[i] = i-1;
        for(int j=1; j <= p[0] && i*p[j] <= n; j++) {
            notp[ i*p[j] ] = 1;
            if(i % p[j] == 0) {phi[ i*p[j] ] = phi[i] * p[j]; break;}
            phi[ i*p[j] ] = phi[i] * (p[j] - 1);
        }
        mod(phi[i] += phi[i-1]);
    }
}

namespace ha {
    const int p = 1001001;
    struct meow{int ne; ll val, r;} e[3000];
    int cnt=1, h[p];
    inline void insert(ll x, ll val) {
        ll u = x % p;
        for(int i=h[u];i;i=e[i].ne) if(e[i].r == x) return;
        e[++cnt] = (meow){h[u], val, x}; h[u] = cnt;
    }
    inline ll quer(ll x) {
        ll u = x % p;
        for(int i=h[u];i;i=e[i].ne) if(e[i].r == x) return e[i].val;
        return -1;
    }
} using ha::insert; using ha::quer;

inline ll sum(ll n) {return n %mo * ((n+1) %mo) %mo *inv2 %mo;}
ll dj_s(ll n) {
    if(n <= U) return phi[n];
    if(quer(n) != -1) return quer(n);
    ll ans = sum(n), r;
    for(ll i=2; i<=n; i=r+1) {
        r = n/(n/i);
        mod(ans -= (r-i+1) %mo * dj_s(n/i) %mo);
    }
    insert(n, ans);
    return ans;
}

ll solve(ll n) {
    ll ans = 0, r;
    for(ll i=1; i<=n; i=r+1) {
        r = n/(n/i);
        mod(ans += dj_s(n/i) * (sum(r) - sum(i-1)) %mo);
    }
    return ans;
}
int main() {
    freopen("in", "r", stdin);
    sieve(U);
    n=read();
    ll ans = 2 * solve(n) %mo - sum(n); mod(ans);
    printf("%lld", ans);
}