题意
给定 \(n\) ,求 \(\sum_{i=1}^n \sum_{j=1}^n lcm(i,j)\)。
\(n\leq 10^{10}\)
分析
- 推式子
\[\begin{aligned} ans &= 2\sum_{i=1}^n\sum_{j=1}^ilcm(i,j)-\sum_{i=1}^nlcm(i,i) \\ &=2\sum_{i=1}^n \sum_{j=1}^i \frac{ij}{(i,j)}-\frac{n(n+1)}{2} \\ &=2\sum_{d=1}^n\frac{1}{d}\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{i}[i\perp j]ijd^2-\frac{n(n+1)}{2} \\&=2\sum_{d=1}^nd\sum_{i=1}^{\frac{n}{d}}i\sum_{j=1}^{i}[i\perp j]j-\frac{n(n+1)}{2}\end{aligned}
\]
又因为 \(\sum_{i=1}^n[i \perp n]i \ =\frac{n\varphi (n)+[n=1]}{2}\)，所以
\[\begin{aligned} ans &=2\sum_{d=1}^nd\sum_{i=1}^{\frac{n}{d}}\frac{i^2\varphi (i)+[i=1]}{2}-\frac{n(n+1)}{2}\\ &=\sum_{d=1}^nd\sum_{i=1}^{\frac{n}{d}}i^2\varphi (i) \end{aligned}
\]
\]
我们记 \(f(i)=i^2\varphi (i), S(n)=\sum_{i=1}^nf(i)\)，那么原式变形为：
\[\sum_{d=1}^nS(\frac{n}{d})d
\]
\]
考虑用杜教筛求 \(S\) ,根据：
\[g(1)S(n)=\sum_{i=1}^n\sum_{j|i}f(j)g(\frac{i}{j})+\sum_{i=2}^ng(i)S(\frac{n}{i})
\]
\]
为了消去 \(i^2\) 的影响，我们考虑构造 \(g(i)=i^2\),得到：
\[\begin{aligned}S(n)&=\sum_{i=1}^n\sum_{j|i}j^2\varphi(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\n", (solve(n) % mod + mod) % mod);
return 0;
}