[51Nod1238]最小公倍数之和 V3[杜教筛]

题意

给定 \(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;
}