【杜教筛】这是一道简单的数学题

题目大意：

求

\[\sum\limits_{i=1}^n\sum\limits^i_{j=1}\frac{\operatorname{lcm}(i,j)}{\gcd(i,j)} \]

\(1\le n\le 10^9\)。

题解：

为了方便，考虑求 \(\sum\limits_{i=1}^n\sum\limits^n_{j=1}\frac{\operatorname{lcm}(i,j)}{\gcd(i,j)}\)

\[\sum\limits_{i=1}^n\sum\limits^n_{j=1}\frac{\operatorname{lcm}(i,j)}{\gcd(i,j)}\\ =\sum\limits_{i=1}^n\sum\limits^n_{j=1}\frac{ij}{\gcd^2(i,j)}\\ =\sum\limits_{d=1}\frac{\sum\limits^n_{i=1}\sum\limits^n_{j=1}ij[\gcd(i,j)=d]}{d^2}\\ =\sum\limits_{d=1}\frac{\sum\limits_{k=1}\mu(k)(\frac{n}{kd})^2(\frac{n}{kd}+1)^2(kd)^2/4}{d^2}\\ =\sum\limits_{d=1}\sum\limits_{k=1}\mu(k)(\frac{n}{kd})^2(\frac{n}{kd}+1)^2k^2/4\\ \]

以下忽略那个 \(/4\)，因为这和整个式子毫无关联。

考虑枚举 \(kd\)，得到原式\(=\)

\[\sum\limits_{d=1}(\frac{n}{d})^2(\frac{n}{d}+1)^2\sum_{k\vert d}\frac{\mu(k)k^2}{d} \]

即

\[\sum\limits_{d=1}(\frac{n}{d})^2(\frac{n}{d}+1)^2\frac{\sum_{k\vert d}\mu(k)k^2}{d} \]

最后一步，考虑使用杜教筛快速求出函数 \(f(n)=\sum_{k\vert d}\mu(k)k^2\) 的前缀和。

令 \(h(n)=\mu(n)n^2\)，那么 \(f(n)=h*1\)。通过观察，发现若构造函数 \(g(n)=n^2\)，则有 \(h*g=\epsilon\)。那么 \(f*g=h*g*1=\epsilon*1=1\)。

综上，我们可以快速求出函数 \(g\) 与函数 \(f*g\) 的前缀和，套杜教筛即可。

#include <cstdio>
#include <cmath>

const int mod = 1e9 + 7, inv2 = 5e8 + 4, inv6 = 166666668;
int f[1500005], Prime[150005], cnt, n;
bool mark[1500005];
class hash_map {
	static const int mod = 19260817;
	int key[10000005], val[10000005], head[10000005], nxt[10000005], tot;
public:
	inline void insert(int x, int y) {
		key[++ tot] = x, val[tot] = y, nxt[tot] = head[x % mod], head[x % mod] = tot;
	}
	inline int find(int x) const {
		for (int i = head[x % mod]; i; i = nxt[i])
			if (key[i] == x) return val[i];
		return -2000000000;
	}
} Sum;

void init(int n) {
	f[1] = 1;
	Sum.insert(0, 0), Sum.insert(1, 1);
	for (int i = 2, sum = 1; i <= n; ++ i) {
		if (!mark[i]) f[i] = (1ll - 1ll * i * i % mod) % mod, Prime[++ cnt] = i;
		Sum.insert(i, sum = (sum + f[i]) % mod);
		for (int j = 1; j <= cnt && i * Prime[j] <= n; ++ j) {
			mark[i * Prime[j]] = true;
			if (i % Prime[j] == 0) {f[i * Prime[j]] = f[i]; break;}
			f[i * Prime[j]] = 1ll * f[i] * f[Prime[j]] % mod;
		}
	}
}

inline long long sq(int n) {return 1ll * n * n % mod * (n + 1) % mod * (n + 1) % mod;}
inline int g(int n) {return 1ll * n * (n + 1) % mod * (2 * n + 1) % mod * inv6 % mod;}
int getsum(int n) {
	if (Sum.find(n) != -2000000000) return Sum.find(n);
	int ans = n;
	for (int l = 2, r; l <= n; l = r + 1) {
		r = n / (n / l);
		ans = (ans - 1ll * (g(r) - g(l - 1)) % mod * getsum(n / l) % mod) % mod;
	}
	Sum.insert(n, ans);
	return ans;
}

int main() {
	init(1500000);
	int n, ans = 0;
	scanf("%d", &n);
	for (int l = 1, r; l <= n; l = r + 1) {	
		r = n / (n / l);
		ans = (ans + sq(n / l) % mod * (getsum(r) - getsum(l - 1)) % mod) % mod;
	}
	printf("%d", ((int)((1ll * ans * inv2 % mod * inv2 % mod + n) % mod * inv2 % mod) + mod) % mod);
	return 0;
}