题意
给定 (n) ,求 (sum_{i=1}^n sum_{j=1}^n lcm(i,j))。
(nleq 10^{10})
分析
- 推式子[egin{aligned} ans &= 2sum_{i=1}^nsum_{j=1}^ilcm(i,j)-sum_{i=1}^nlcm(i,i) \ &=2sum_{i=1}^n sum_{j=1}^i frac{ij}{(i,j)}-frac{n(n+1)}{2} \ &=2sum_{d=1}^nfrac{1}{d}sum_{i=1}^{frac{n}{d}}sum_{j=1}^{i}[iperp j]ijd^2-frac{n(n+1)}{2} \&=2sum_{d=1}^ndsum_{i=1}^{frac{n}{d}}isum_{j=1}^{i}[iperp j]j-frac{n(n+1)}{2}end{aligned} ]
又因为 (sum_{i=1}^n[i perp n]i =frac{nvarphi (n)+[n=1]}{2})，所以
[egin{aligned} ans &=2sum_{d=1}^ndsum_{i=1}^{frac{n}{d}}frac{i^2varphi (i)+[i=1]}{2}-frac{n(n+1)}{2}\ &=sum_{d=1}^ndsum_{i=1}^{frac{n}{d}}i^2varphi (i) end{aligned}
]
我们记 (f(i)=i^2varphi (i), S(n)=sum_{i=1}^nf(i))，那么原式变形为：
[sum_{d=1}^nS(frac{n}{d})d
]
考虑用杜教筛求 (S) ,根据：
[g(1)S(n)=sum_{i=1}^nsum_{j|i}f(j)g(frac{i}{j})+sum_{i=2}^ng(i)S(frac{n}{i})
]
为了消去 (i^2) 的影响，我们考虑构造 (g(i)=i^2),得到：
[egin{aligned}S(n)&=sum_{i=1}^nsum_{j|i}j^2varphi(j)frac{i^2}{j^2}+sum_{i=2}^ni^2S(frac{n}{i}) \ &=sum_{i=1}^ni^3+sum_{i=2}^ni^2S(frac{n}{i})end{aligned}
]
杜教筛上式，外层数论分块即可。
复杂度就不算了哈
代码
#include<bits/stdc++.h>
using namespace std;
typedef long long LL;
#define go(u) for(int i = head[u], v = e[i].to; i; i=e[i].lst, v=e[i].to)
#define rep(i, a, b) for(int i = a; i <= b; ++i)
#define pb push_back
inline int gi() {
int x = 0,f = 1;
char ch = getchar();
while(!isdigit(ch)) {
if(ch == '-') f = -1;
ch = getchar();
}
while(isdigit(ch)) {
x = (x << 3) + (x << 1) + ch - 48;
ch = getchar();
}
return x * f;
}
template <typename T> inline void Max(T &a, T b){if(a < b) a = b;}
template <typename T> inline void Min(T &a, T b){if(a > b) a = b;}
const LL mod = 1e9 + 7, N = 2e6 + 7;
LL n, pc, inv2 = 500000004, inv6;
LL f[N], pri[N];
bool vis[N];
map<LL, LL> F;
LL Pow(LL a, LL b) {
LL res = 1ll;
for(; b; b >>= 1, a = a * a %mod) if(b & 1) res = res * a % mod;
return res;
}
void pre(int n) {
f[1] = 1;int to;
for(int i = 2; i <= n; ++i) {
if(!vis[i]) { pri[++pc] = i, f[i] = 1ll * i * i % mod * (i - 1) % mod; }
for(int j = 1; j <= pc && (to = i * pri[j]) <= n; ++j) {
vis[to] = 1;
if(i % pri[j] == 0) {
f[to] = f[i] * pri[j] % mod * pri[j] % mod * pri[j] %mod;
break;
}
f[to] = f[i] * f[pri[j]] % mod;
}
}
for(int i = 2; i <= n; ++ i) f[i] = (f[i - 1] + f[i]) % mod;
}
LL s3(LL n) {
n %= mod;
LL x = (n + 1) * n % mod * inv2 %mod;
return x * x % mod;
}
LL s2(LL n) {
n %= mod;
return n * (n + 1) % mod * (2 * n + 1) % mod * inv6 % mod;
}
LL s1(LL n) {
n %= mod;
return n * (n + 1) % mod * inv2 % mod;
}
void add(LL &a, LL b){ a += b; if(a >= mod) a -= mod;}
LL gg(LL n) {
if(n <= 2e6) return f[n];
if(F.count(n)) return F[n];
LL res = s3(n);
for(LL i = 2, lst = 2; i <= n; i = lst + 1) {
lst = n / (n / i);
add(res, mod - (s2(lst) - s2(i - 1) + mod) * gg(n / i) %mod);
}
return F[n] = res;
}
LL solve(LL n) {
LL res = 0ll;
for(LL i = 1, lst = 1; i <= n; i = lst + 1) {
lst = n / (n / i);
add(res, (s1(lst) - s1(i - 1) + mod) * gg(n / i) % mod);
}
return res;
}
int main() {
scanf("%lld", &n);
inv6 = Pow(6, mod - 2);
pre(2e6);
printf("%lld
", (solve(n) % mod + mod) % mod);
return 0;
}