题解
参考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;
}