51nod 1237 杜教筛

题解

参考1238

code:

#include <bits/stdc++.h>using namespace std;typedef long long ll;const int N = 1000001;const ll mod = 1e9 + 7;const ll inv = (mod + 1) / 2;const int mo = 2333333;bool isPrime[N];ll phi[N];int prime[N];int cnt;void init(){    memset(isPrime, true, sizeof isPrime);    phi[1] = 1;    for(int i = 2; i < N; ++i)    {        if(isPrime[i])        {            prime[++cnt] = i;            phi[i] = i - 1;        }        for(int j = 1; j <= cnt && i * prime[j] < N; ++j)        {            isPrime[i * prime[j]] = false;            if(i % prime[j] == 0)            {                phi[i * prime[j]] = phi[i] * prime[j];                break;            }            phi[i * prime[j]] = phi[i] * (prime[j] - 1);        }    }    for(int i = 2; i < N; ++i)        phi[i] = (phi[i] + phi[i - 1]) % mod;}ll n;int last[mo], next[mo];ll t[mo], v[mo];int l;void add(int x, ll y, ll z){    t[++l] = y;    next[l] = last[x];    last[x] = l;    v[l] = z;}ll cal(ll x){    if(x < N) return phi[x];    int k = x % mo;    for(int i = last[k]; i; i = next[i])        if(t[i] == x) return v[i];    ll res = x % mod * ((x + 1) % mod) % mod * inv % mod;    ll r;    for(ll i = 2; i <= x; i = r + 1)    {        ll tmp = x / i;        r = x / tmp;        res = ((res - ((r - i + 1) % mod * cal(tmp) % mod)) % mod + mod) % mod;    }    add(k, x, res);    return res;}ll Calc(ll x){    ll res = 0;    ll r;    for(ll i = 1; i <= x; i = r + 1)    {        ll tmp = x / i;        r = x / tmp;        res = (res + ((r - i + 1) % mod * ((r + i) % mod) % mod * inv % mod * cal(tmp) % mod)) % mod;    //    cout << res << endl;    }    res = res * 2 % mod;    ll t = x % mod * ((x + 1) % mod) % mod * inv % mod;    return ((res - t) % mod + mod) % mod;}int main(){    init();    cin >> n;    cout << Calc(n) << endl;    return 0;}