传送门:51nod 1238
题意
求
G(N)=∑i=1N∑j=1Nlcm(i,j)
题解
首先
G(N)=∑i=1N∑j=1Nlcm(i,j)=2∑i=1N∑j=1ilcm(i,j)−∑i=1Nlcm(i,i)=2∑i=1Ni∑d|i∑u=1idu[gcd(u,id)=1]−N(N+1)2=2∑i=1Ni∑d|iidϕ(id)+[id=1]2−N(N+1)2=∑i=1Ni∑d|iidϕ(id)+[id=1]−N(N+1)2=∑i=1Ni(1+∑d|iidϕ(id))−N(N+1)2=∑i=1Ni∑d|iidϕ(id)=∑i=1Ni∑d|idϕ(d)
其次枚举因子对i的贡献:
G(N)=∑i=1Ni∑d|idϕ(d)=∑i=1N∑d=1⌊Ni⌋dϕ(d)
可用分块:
令
F(x)=∑i=1xiϕ(i)
预处理 N23 以内的 F(i)
对于 x>N23
杜教筛求解:
F(x)=∑i=1xiϕ(i)=∑i=1xi(i−∑d|i,d<iϕ(d))=∑i=1xi2−∑i=1xi∑d|i,d<iϕ(d))=x(2x+1)(x+1)6−∑i=2xi∑d=1⌊xi⌋dϕ(d)=x(2x+1)(x+1)6−∑i=2xiF(⌊xi⌋)
分块求解:
code:
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const ll mod = 1e9 + 7;
const ll inv = (mod + 1) / 2;
const ll inv2 = inv * (mod + 1) / 3 % mod;
const int mo = 2333333;
const int N = 1000001;
bool isPrime[N];
int prime[N];
ll phi[N];
ll sum[N];
void init()
{
int cnt = 0;
memset(isPrime, true, sizeof isPrime);
isPrime[1] = false;
phi[1] = 1;
sum[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(ll i = 2; i < N; ++i)
{
sum[i] = (sum[i - 1] + i * i * phi[i] % mod) % mod;
}
}
ll Sum(ll x)
{
return x % mod * ((x * 2 + 1) % mod) % mod * ((x + 1) % mod) % mod * inv2 % mod;
}
ll getSum(ll R, ll L)
{
return ((Sum(R) - Sum(L)) % mod + mod) % mod;
}
ll n;
int next[mo], last[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 sum[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;
res = res * res % mod;
// if(x == 10000000000) cout << res << endl;
ll r;
for(ll i = 2; i <= x; i = r + 1)
{
ll tmp = x / i;
r = x / tmp;
ll mul = getSum(r, i - 1);
//if(x == 10000000000)cout << mul << " ";
ll ans = cal(tmp) * mul % mod;
// if(x == 10000000000)cout << ans << " ";
res = ((res - ans) % mod + mod) % mod;
// if(x == 10000000000) cout << res << endl;
}
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;
ll mul = ((r - i + 1) % mod) % mod * ((r + i) % mod) % mod * inv % mod;
// cout << mul << " ";
res = (res + mul * cal(tmp) % mod) % mod;
//cout << res << endl;
}
r = 0;
/** r = x * inv % mod;
ll mul = 1;
for(ll i = 1; i * i <= x; ++i)
{
if(x % i == 0)
{
mul = (mul + i * getPhi(i)) % mod;
if(i * i != x) mul = (mul + x / i * getPhi(x / i)) % mod;
}
}
r = r * mul % mod;*/
return ((res - r) % mod + mod) % mod;
}
int main()
{
init();
cin >> n;
// cout << inv << " " <<inv2 << endl;
cout << Calc(n) << endl;
return 0;
}