筛法

杜教筛，$\textsf {Min-25}$ 筛

越来越边缘了啊，主要是板子，题很少

建议阅读 数论函数基础，数论分块，莫比乌斯函数，欧拉函数 章节，了解基本概念与先要知识

全文 绝大多数 内容是对 [0] 中讲述的 粗略抄写 和 胡乱加工

1. 杜教筛

在 [0] 中被 跳过了，悲，所以这部分主要参考 [1] $\color {black} \textsf {harryzhr}$ 老师 $2023$ 年的一份课件

但是课件中 时间复杂度 部分是 伪证一例，故此处参考 [2] OI - Wiki 上的证明

概念

对于一个 数论函数 $f$，杜教筛 可以在 低于线性 的时间复杂度内计算 $S (n) = \sum _ {i = 1} ^ n f (i)$

也就是 更快的计算求和函数

做法

考虑找一个 积性函数 $g$，我们考虑 $f, g$ 狄利克雷卷积 的 求和函数，即

\[\begin {aligned} \sum _ {i = 1} ^ n (f * g) (i) \end {aligned} \]

简单转化，有

\[\begin {aligned} \sum _ {i = 1} ^ n (f * g) (i) &= \sum _ {i = 1} ^ n { \sum _ {d \mid i} { g (d) f \left ( \dfrac i d \right ) } } \\ &= \sum _ {i = 1} ^ n { g (i) S \left ( \left \lfloor \dfrac n i \right \rfloor \right ) } \end {aligned} \]

由于 $g$ 为 积性函数，故有 $g (1) = 1$，于是 $S (n) = g(1)S(n)$，有下式

\[\begin {aligned} g(1)S(n) &= \sum _ {i = 1} ^ n { g (i) S \left ( \left \lfloor \dfrac n i \right \rfloor \right ) } - \sum _ {i = 2} ^ n { g (i) S \left ( \left \lfloor \dfrac n i \right \rfloor \right ) } \\ &= \sum _ {i = 1} ^ n { (f * g) (i) } - \sum _ {i = 2} ^ n { g (i) S \left ( \left \lfloor \dfrac n i \right \rfloor \right ) } \end {aligned} \]

设 $f * g = h$，那么如果 $h, g$ 的 前缀和 是容易求的，那么使用 数论分块 就可递归求解 $\sum f$

这个东西简单来讲其实就是直到 $g, f * g$ 的 求和函数，快速求解 $f$ 的 求和函数

时间复杂度

若 $h, g$ 的 前缀和 均可 $O (1)$ 取得，由 数论分块 性质，我们知道 $\left \lfloor \dfrac n i \right \rfloor$ 有 $O (\sqrt n)$ 种取值

我们不妨设取值集合为 $R (n) = \left \{ \left \lfloor \dfrac n k \right \rfloor : k \in [1, n] \right \}$，存在如下引理

对于任意的 $m \in R(n)$，有 $R (m) \subseteq R (n)$

证明是容易的，考虑若 $m = \left \lfloor \dfrac n x \right \rfloor$，则任意 $\left \lfloor \dfrac m y \right \rfloor = \left \lfloor \dfrac n {xy} \right \rfloor \in R (n)$

这告诉我们，所有递归过程中涉及到的 $S (x)$ 中的 $x$，取值均被包含在 $R (n)$ 中

而我们若 使用记忆化，则对于每个 $S (x)$，我们至多计算一次

换言之，我们至多对于 $R (n)$ 中每个值 $x$ 计算 $S (x)$ 就能得到最终答案，于是

\[\begin {aligned} T (n) &= \sum _ {k \in R (n)} T (k) \\ &= O (\sqrt n) + \sum _ {k = 1} ^ { \left \lfloor \sqrt n \right \rfloor } O (\sqrt k) + \sum _ {k = 2} ^ { \left \lfloor \sqrt n \right \rfloor } O \left ( \sqrt { \dfrac n k } \right ) \\ &= O \left ( \int _ 0 ^ {\sqrt n} { \left ( \sqrt x + \sqrt {\dfrac n x} \right ) dx } \right ) \\ &= O (n ^ {0.75}) \end {aligned} \]

而如果我们能在合理的时间 $(T _ 0 (m))$ 内，预处理出一部分 $S (k) ~ (k \in [1, m], m \ge \left \lfloor \sqrt n \right \rfloor)$

则时间复杂度将变成

\[\begin {aligned} T (n) &= T _ 0 (m) + \sum _ {k \in R (n), k > m} T (k) \\ &= T _ 0 (m) + \sum _ {k = 1} ^ { \left \lfloor \frac n m \right \rfloor } O \left ( \sqrt { \dfrac n k } \right ) \\ &= O \left ( T_0 (m) + \int _ {0} ^ {\frac n m} { \sqrt {\dfrac n x} dx } \right ) \\ &= O \left ( T_0 (m) + \dfrac n {\sqrt m} \right ) \end {aligned} \]

若使用线性筛，即 $T _ 0 (m) = O (m)$，则由 均值不等式，当 $m = n ^ {\frac 2 3}$ 时， $T (n) _ {\min} = O (n ^ \frac 2 3)$

Luogu P4213 【模板】杜教筛

求 $S_1 (n) = \sum _ {i = 1} ^ n \varphi (i)$，$S_2 (n) = \sum _ {i = 1} ^ n \mu (i)$

对于 $S_1$，考虑有 $\varphi * 1 = \operatorname {id}$，故有

\[\begin {aligned} S_1 (n) &= \sum _ {i = 1} ^ n { \operatorname {id} (i) } - \sum _ {i = 2} ^ n { S_1 \left ( \left \lfloor \dfrac n i \right \rfloor \right ) } \\ &= \dfrac {n (n + 1)} 2 - \sum _ {i = 2} ^ n { S_1 \left ( \left \lfloor \dfrac n i \right \rfloor \right ) } \end {aligned} \]

对于 $S_2$，考虑有 $\mu * 1 = \epsilon$，故有

\[\begin {aligned} S_2 (n) &= \sum _ {i = 1} ^ n { \epsilon (i) } - \sum _ {i = 2} ^ n { S_2 \left ( \left \lfloor \dfrac n i \right \rfloor \right ) } \\ &= 1 - \sum _ {i = 2} ^ n { S_2 \left ( \left \lfloor \dfrac n i \right \rfloor \right ) } \end {aligned} \]

然后递归求解即可，记得需要 记忆化

悲惨的消息是 std::unordered_map 常数太大，所以这个题 $O (n ^ {0.75})$ 的方法是 过不去的

于是我们只能 线性筛预处理部分前缀和，然后做到 $O (n ^ {\frac 2 3})$ 的复杂度，才能通过

#include <bits/stdc++.h>

const int MAXN = 1500005;
const int maxn = 1500000;

std::unordered_map <int32_t, int64_t> sum1, sum2;
int32_t vis[MAXN], pri[MAXN], phi[MAXN], mu[MAXN], cnt;
int64_t s1[MAXN], s2[MAXN];

inline void Sieve (const int n = maxn) {
	phi[1] = mu[1] = 1;
	for (int i = 2; i <= n; ++ i) {
		if (!vis[i]) pri[++ cnt] = i, mu[i] = -1, phi[i] = i - 1;
		for (int k = 1; k <= cnt && pri[k] * i <= n; ++ k) {
			vis[i * pri[k]] = 1;
			if (i % pri[k] == 0) {
				phi[i * pri[k]] = pri[k] * phi[i];
				break ;
			}
			mu[i * pri[k]] = - mu[i];
			phi[i * pri[k]] = (pri[k] - 1) * phi[i]; 
		}
	}
	for (int i = 1; i <= n; ++ i) s1[i] = s1[i - 1] + phi[i];
	for (int i = 1; i <= n; ++ i) s2[i] = s2[i - 1] + mu [i];
}

inline int64_t Sieve1 (const int64_t n) {
	if (n <= maxn) return s1[n];
	if (sum1.count (n)) return sum1[n];
	int64_t ans = n * (n + 1) / 2;
	for (int64_t l = 2, r; l <= n; l = r + 1) 
		r = n / (n / l), ans = ans - (r - l + 1) * Sieve1 (n / l);
	return sum1[n] = ans;
}

inline int64_t Sieve2 (const int64_t n) {
	if (n <= maxn) return s2[n];
	if (sum2.count (n)) return sum2[n];
	int64_t ans = 1;
	for (int64_t l = 2, r; l <= n; l = r + 1) 
		r = n / (n / l), ans = ans - (r - l + 1) * Sieve2 (n / l);
	return sum2[n] = ans;
}

int32_t N;
int64_t x;

int main () {
	
	Sieve (), sum1[0] = sum2[0] = 0;
	
	std::cin >> N;
	
	for (int i = 1; i <= N; ++ i) {
		std::cin >> x;
		std::cout << Sieve1 (x) << ' ' << Sieve2 (x) << '\n';
	}
	
	return 0;
}

Luogu P3768 简单的数学题

并不困难，我们考虑先对式子进行转化，惯例枚举 $\gcd$

\[\begin {aligned} \sum _ {d = 1} ^ n { \sum _ {i = 1} ^ n { \sum _ {j = 1} ^ n { [\gcd (i, j) = d] ij } } } \end {aligned} \]

只枚举 $d$ 的倍数

\[\begin {aligned} \sum _ {d = 1} ^ n { d \sum _ {i = 1} ^ { \left \lfloor \frac n d \right \rfloor } { \sum _ {j = 1} ^ { \left \lfloor \frac n d \right \rfloor } { [\gcd (i, j) = 1] ijd ^ 2 } } } \end {aligned} \]

$d ^ 2$ 放到前面去，然后使用 莫比乌斯函数 - Trick 3 - ex，有

\[\begin {aligned} \sum _ {d = 1} ^ n { d ^ 3 \sum _ {i = 1} ^ { \left \lfloor \frac n d \right \rfloor } i ^ 2 \varphi (i) } \end {aligned} \]

这时候已经结束了，考虑如果我们能快速求得 后半部分，那么整除分块就可以做

于是考虑 杜教筛，筛出 $S (n) = \sum _ {i = 1} ^ n i ^ 2 \varphi (i)$ 的值，注意到 欧拉反演 $\sum _ {d \mid n} \varphi (d) = n$

我们有 $f (i) = i ^ 2 \varphi (i)$，可以想到设 $g (i) = i ^ 2$，于是

\[\begin {aligned} h (n) &= (f * g) (n) \\ &= \sum _ {d \mid n} { f (d) g \left ( \dfrac n d \right ) } \\ &= \sum _ {d \mid n} { d ^ 2 \varphi (d) \left ( \dfrac n d \right ) ^ 2 } \\ &= n ^ 3 \end {aligned} \]

牛得很，于是就可以做了

\[\begin {aligned} S (n) &= \sum _ {i = 1} ^ n h (i) - \sum _ {i = 2} ^ n { g (i) S \left ( \left \lfloor \dfrac n i \right \rfloor \right ) } \\ &= \sum _ {i = 1} ^ n i ^ 3 - \sum _ {i = 2} ^ n { i ^ 2 S \left ( \left \lfloor \dfrac n i \right \rfloor \right ) } \end {aligned} \]

需求 $\sum i ^ 2$ 与 $\sum i ^ 3$ 的公式

但是杜教筛一次不是只能求 求和函数单点值 吗？

如果这样那么时间复杂度就是 $O (n ^ {0.25} n ^ {\frac 2 3})$，爆掉了

实质上，筛一次我们求出了 $R (n)$ 这个集合内所有点的值，而这正是我们这个题需要的

故时间复杂度 $O (n ^ {\frac 2 3} + \sqrt n)$，可以通过本题

#include <bits/stdc++.h>

const int MAXN = 5000005;
const int PREN = 5000000;

int64_t MOD = 998244353;

std::unordered_map <int64_t, int64_t> sum;

int64_t s[MAXN], I2, I6;
int64_t vis[MAXN], pri[MAXN], phi[MAXN], cnt;

inline int64_t qpow (int64_t a, int64_t b = MOD - 2) {
	int64_t res = 1;
	a = a % MOD;
	while (b) {
		if (b & 1) res = res * a % MOD;
		a = a * a % MOD, b >>= 1;
	}
	return res;
}

inline void Sieve (const int64_t n = PREN) {
	phi[1] = 1, I2 = qpow (2), I6 = qpow (6);
	for (int64_t i = 2; i <= n; ++ i) {
		if (!vis[i]) pri[++ cnt] = i, phi[i] = i - 1;
		for (int64_t k = 1; k <= cnt && pri[k] * i <= n; ++ k) {
			vis[i * pri[k]] = 1;
			if (i % pri[k] == 0) {
				phi[i * pri[k]] = pri[k] * phi[i];
				break ;
			}
			phi[i * pri[k]] = (pri[k] - 1) * phi[i]; 
		}
	}
	for (int64_t i = 1; i <= n; ++ i) s[i] = 1ll * i * i % MOD * phi[i] % MOD;
	for (int64_t i = 1; i <= n; ++ i) s[i] = (s[i - 1] + s[i]) % MOD;
} 

inline int64_t Sum1 (int64_t x) {
	x = x % MOD;
	return x * (x + 1) % MOD * I2 % MOD;
}

inline int64_t Sum2 (int64_t x) {
	x = x % MOD;
	return x % MOD * (x + 1) % MOD * (x + x + 1) % MOD * I6 % MOD;
}

inline int64_t Sum3 (const int64_t x) {
	return Sum1 (x) * Sum1 (x) % MOD;
}

inline int64_t Sieve1 (const int64_t n) {
	if (n <= PREN) return s[n];
	if (sum.count (n)) return sum[n];
	int64_t ans = Sum3 (n);
	for (int64_t l = 2, r; l <= n; l = r + 1) 
		r = n / (n / l), ans = (ans - (Sum2 (r) - Sum2 (l - 1) + MOD) % MOD * Sieve1 (n / l) % MOD + MOD) % MOD; 
	return sum[n] = ans;
}

inline int64_t Solve (const int64_t n) {
	int64_t ans = 0;
	for (int64_t l = 1, r; l <= n; l = r + 1) 
		r = n / (n / l), ans = (ans + (((Sum3 (r) - Sum3 (l - 1)) % MOD + MOD) % MOD * Sieve1 (n / l) % MOD + MOD) % MOD) % MOD;
	return ans;
}

int64_t n;

int main () {
	
	std::cin >> MOD >> n;
	
	Sieve (), Sieve1 (n);
	
	std::cout << Solve (n) << '\n';
	
	return 0;
}

Luogu P1587 [NOI2016] 循环之美

题目大意就是计数在 $k$ 进制下 纯循环小数 的取值个数

我们来分析题目中给出的条件，首先 每个取值只计算一次，我们可以转化成 计数最简分数

即对于分数 $\dfrac x y$，要求 $x \perp y$

然后考虑刻画 纯循环小数 这个性质，注意到其满足 小数部分左 / 右移 $l$ 位后不变

即 $x \equiv x k ^ l \pmod y$（其中 $k ^ l$ 表示左移 $l$ 位，$\bmod y$ 表示取小数部分）

于是有 $k ^ l \equiv 1 \pmod y$，我们可以证明这个条件等价于 $k \perp y$

设 $\gcd (k, y) = d$，由于 $k ^ l \equiv 1 \pmod y$，故有 $k ^ l = my + 1$（$m \in \mathbb Z$）

由于 $d \mid k \wedge d \mid y$，显然需要有 $d \mid my + 1$，即 $d \mid 1$，故 $d = 1$，$k \perp y$

所以这个题就等价于让我们计数满足 $x \perp y$ 且 $k \perp y$ 的 二元组个数，即下式

\[\begin {aligned} \sum _ {x = 1} ^ n { \sum _ {y = 1} ^ m { [\gcd (x, y) = 1] [\gcd (y, k) = 1] } } \end {aligned} \]

容易利用 莫比乌斯反演，转换枚举顺序 后得到

\[\begin {aligned} \sum _ {d = 1} ^ {\min (n, m)} { \sum _ {x = 1} ^ { \left \lfloor \frac n d \right \rfloor } { \sum _ {y = 1} ^ { \left \lfloor \frac m d \right \rfloor } { \mu (d) [\gcd (yd, k) = 1] } } } \end {aligned} \]

将无用的 $x$ 提出，分离 $y, d$，于是

\[\sum _ {d = 1} ^ {\min (n, m)} { \mu (d) \left \lfloor \dfrac n d \right \rfloor [\gcd (d, k) = 1] } \sum _ {y = 1} ^ { \left \lfloor \frac m d \right \rfloor } { [\gcd (y, k) = 1] } \]

前面部分的处理没有头绪，但是后面一个 $\sum$ 部分是简单的，设 $g (n) = \sum _ {y = 1} ^ n [\gcd (i, k) = 1]$

显然这东西 $[\gcd (i, k) = 1]$ 这东西有循环节长 $k$，容易知道其通项公式

\[\begin {aligned} g (n) &= \left \lfloor \dfrac n k \right \rfloor \varphi (k) + g (n \bmod k) \\ &= \left \lfloor \dfrac n k \right \rfloor g (k) + g (n \bmod k) \end {aligned} \]

于是我们只需要 $O (k)$ 预处理 $1 \sim k$ 的答案即可 $O (1)$ 求解

然后此时我们可以将要求的式子转化成下面的形式

\[\begin {aligned} \sum _ {d = 1} ^ {\min (n, m)} { \left \lfloor \dfrac n d \right \rfloor g \left ( \left \lfloor \dfrac m d \right \rfloor \right ) \mu (d) [\gcd (d, k) = 1] } \end {aligned} \]

显然前两项可以 整除分块 处理掉，于是我们只在意 $f (n) = \sum _ {i = 1} ^ n \mu (i) [\gcd (i, k) = 1]$ 这部分

用类似 杜教筛 的式子，我们可以把它表示成 两和式相减 的形式

\[\begin {aligned} f (n) = \sum _ {i = 1} ^ n { [\gcd (i, k) = 1] f \left ( \left \lfloor \dfrac n i \right \rfloor \right ) } - \sum _ {i = 2} ^ n { [\gcd (i, k) = 1] f \left ( \left \lfloor \dfrac n i \right \rfloor \right ) } \end {aligned} \]

前面部分展开就变成了（下式不包含后面部分）

\[\begin {aligned} \sum _ {i = 1} ^ n { [\gcd (i, k) = 1] \sum _ {j = 1} ^ { \left \lfloor \frac n i \right \rfloor } { \mu (j) [\gcd (j, k) = 1] } } \end {aligned} \]

把 $\sum _ {j}$ 放到前面，合并 $\gcd$ 条件，有

\[\begin {aligned} \sum _ {i = 1} ^ n { \sum _ {j = 1} ^ { \left \lfloor \frac n i \right \rfloor } { \mu (j) [\gcd (ij, k) = 1] } } \end {aligned} \]

经典套路枚举 $d = ij$，$j$ 作为因数

\[\begin {aligned} \sum _ {d = 1} ^ n { \sum _ {j \mid d} { \mu (j) [\gcd (d, k) = 1] } } \end {aligned} \]

我们在其中发现了 $\sum _ {j \mid d} \mu (j)$，好得很

\[\begin {aligned} \sum _ {d = 1} ^ n { [d = 1] [\gcd (d, k) = 1] } \end {aligned} \]

于是这部分就是 $1$，那么 $f (n)$ 就可以表示为

\[\begin {aligned} 1 - \sum _ {i = 2} ^ n { [\gcd (i, k) = 1] f \left ( \left \lfloor \dfrac n i \right \rfloor \right ) } \end {aligned} \]

这就是一个 典型的杜教筛 形式了，我们筛出 $f (n)$，套上前面 整除分块 即可

时间复杂度似乎带点玄学，因为 整除分块 致使我们需要 $O (\sqrt n)$ 次 杜教筛

显然 $O (n ^ {\frac 1 2 + \frac 2 3})$ 的复杂度 怎么想都爆了

但显然这个题仍然具有上一个题 筛一次得到多个有效值 的特性

只是在这个题中，具体 跑满筛了多少次 似乎较难刻画，所以说是带点玄学的

#include <bits/stdc++.h>

const int MAXN = 1000005;
const int PERN = 1000000;

int vis[MAXN], pri[MAXN], mu[MAXN], g[MAXN], f[MAXN], cnt;
int N, M, K;

std::unordered_map <int32_t, int64_t> sum;

inline void Sieve (const int n = PERN) {
	for (int i = 1; i <= n; ++ i) g[i] = g[i - 1] + (std::__gcd (i, K) == 1);
	mu[1] = f[1] = 1;
	for (int i = 2; i <= n; ++ i) {
		if (!vis[i]) pri[++ cnt] = i, mu[i] = - 1;
		for (int k = 1; k <= cnt && pri[k] * i <= n; ++ k) {
			vis[i * pri[k]] = 1;
			if (i % pri[k] == 0) break ;
			mu[pri[k] * i] = - mu[i];
		}
		f[i] = f[i - 1] + mu[i] * (g[i] - g[i - 1]);
	}
}

inline int64_t G (const int x) {
	return x / K * g[K] + g[x % K];
}

inline int64_t F (const int n) {
	if (n <= PERN) return f[n];
	if (sum.count (n)) return sum[n];
	int64_t ans = 1;
	for (int l = 2, r; l <= n; l = r + 1) 
		r = n / (n / l), ans = ans - F (n / l) * (G (r) - G (l - 1));
	return sum[n] = ans;
}

inline int64_t Solve (const int n, const int m) {
	int64_t ans = 0;
	for (int l = 1, r; l <= std::min (n, m); l = r + 1) 
		r = std::min (n / (n / l), m / (m / l)), ans = ans + (F (r) - F (l - 1)) * (n / l) * G (m / l);
	return ans;
}

int main () {
	
	std::cin >> N >> M >> K, Sieve ();
	
	std::cout << Solve (N, M) << '\n';
	
	return 0;
}

2. $\textsf {Min-25}$ 筛

概念

快速筛 积性函数 $f$ 前缀和的算法，要求 $f (p ^ k) ~ (p \in \mathbb P)$ 为关于 $p ^ k$ 的 常数次多项式

也就是 $f (p ^ k)$ 容易求得

做法

算法流程简单来说就是先 构造完全积性函数 $f'$ 使得 $\forall p ^ k, p \in \mathbb P, k \in \mathbb Z, f (p ^ k) = f' (p ^ k)$

通过函数 $f'$，我们可以推出 $f$ 在所有 质数幂 处的取值

然后构造函数 $S$，利用 质数幂处的取值 以及 $g$ 函数的值求解 原函数前缀和

时间复杂度

有说 $O \left (\dfrac {n ^ {0.75}} {\log n} \right )$ 的，但事实上在 $OI$ 数据范围里可能比 杜教筛 的 $O (n ^ {\frac 2 3})$ 更快一些

常数较小，可以通过 $10 ^ {10}$ 的数据，具体证明笔者并不会，可以参考 [2] OI Wiki 上的分析

Luogu P5325 【模板】Min_25 筛

此处参考 [3]

设 完全积性函数 $f'$ 满足上述条件，我们钦定一个新函数 $g (n, p)$，有

\[\begin {aligned} g (n, j) = \sum _ {i = 1} ^ n f'(i) (i \in \mathbb P \vee \min _ {p \mid i} p > P _ j) \end {aligned} \]

$P_j$ 表示第 $j$ 个质数

那这玩意儿怎么求呢？考虑通过 $g (n, j - 1)$ 来 转移得出

考虑 $g (n, j)$ 与 $g (n, j - 1)$ 的差距

容易发现，$g (n, j - 1)$ 比 $g (n, j)$ 多的就是 最小质因数 为 $P_j$ 的 $i$ 对应的 $f' (i) $ 的值

我们可以 强行钦定 这些数一定有质因子 $P_j$，于是把这些数除以 $P_j$

而这些数剩下部分的 最小质因子 需要 大于等于 $P_j$，于是这些数贡献就是 $g \left (\left \lfloor \dfrac n {P_j} \right \rfloor, j - 1 \right)$

然后我们还需要 去掉小于 $P_j$ 的质数 的贡献，于是需要减去 $\sum _ {i = 1} ^ {j - 1} f' (P _ i)$

\[\begin {aligned} g (n, j) = g (n, j - 1) - f' (P _ j) \left ( g \left ( \dfrac n {P _ j} , j - 1 \right ) - \sum _ {i = 1} ^ {j - 1} f' (P_i) \right ) \end {aligned} \]

考虑求值，设 $S (n, j) = \sum _ {i = 1} ^ n f (i) (\min _ {p \mid i} p > P _ j)$，我们分成质数和合数来计算贡献

在质数时是简单的，考虑 $f'$ 与 $f$ 在 质数幂 处取值相同，故质数时答案就是

\[\begin {aligned} g \left ( n , |P| \right ) - \sum _ {i = 1} ^ j f' (P _ i) \end {aligned} \]

合数时我们可以 枚举每个数的最小质因数以及次数，根据积性函数性质转移，即

\[\begin {aligned} \sum _ {k > j} { \sum _ {c = 1} ^ {P _ k ^ c \le n} { f (P _ k ^ c) \left ( S \left ( \dfrac n {P _ k ^ c} , k \right ) + [c > 1] \right ) } } \end {aligned} \]

这里 $c > 1$ 的意义是显然的，我们 不要重复统计质数的贡献

但是我们需要统计 质数幂（但不是质数，也就是说不是一次）的贡献

于是有 $S$ 的转移式

\[\begin {aligned} S (n, j) = g (n, |P|) - \sum _ {i = 1} ^ j f' (P_i) + \sum _ {k > j} { \sum _ {c = 1} ^ {P _ k ^ c \le n} { f (P _ k ^ c) \left ( S \left ( \dfrac n {P _ k ^ c} , k \right ) + [c > 1] \right ) } } \end {aligned} \]

最终答案就是 $S (n, 0)$，但是加上 $f (1)$，因为 $S$ 函数未统计 $1$ 处的贡献

#include <bits/stdc++.h>

const int MAXN = 200005;
const int MOD  = 1000000007;
const int I2   = 500000004;
const int I6   = 166666668;

int64_t vis[MAXN], pri[MAXN], cnt = 0;
int64_t pos[MAXN], id1[MAXN], id2[MAXN], tot = 0, sqrtn;
int64_t N;
int64_t h1[MAXN], h2[MAXN], g1[MAXN], g2[MAXN];

inline void Sieve (const int n = sqrtn) {
	for (int i = 2; i <= n; ++ i) {
		if (!vis[i]) pri[++ cnt] = i;
		for (int k = 1; k <= cnt && pri[k] * i <= n; ++ k) {
			vis[pri[k] * i] = 1;
			if (i % pri[k] == 0) break ;
		}
	}
	for (int i = 1; i <= cnt; ++ i) h1[i] = (h1[i - 1] + pri[i]) % MOD;
	for (int i = 1; i <= cnt; ++ i) h2[i] = (h2[i - 1] + pri[i] * pri[i] % MOD) % MOD;
}

inline int32_t Id (const int64_t x) {
	return x <= sqrtn ? id1[x] : id2[N / x];
}

inline int64_t S1 (int64_t x) {
	x = x % MOD;
	return x * (x + 1) % MOD * I2 % MOD;
}

inline int64_t S2 (int64_t x) {
	x = x % MOD;
	return x * (x + 1) % MOD * (x + x + 1) % MOD * I6 % MOD;
}

inline int64_t Fp (int64_t n) {
	n = n % MOD;
	return n * (n - 1) % MOD;
}

inline int64_t G (const int64_t n) {
	return (g2[Id (n)] - g1[Id (n)] + MOD) % MOD;
}

inline int64_t H (const int64_t n) {
	return (h2[n] - h1[n] + MOD) % MOD;
}

inline int64_t S (const int64_t n, const int k) {
	if (pri[k] > n) return 0;
	int64_t ans = (G (n) - H (k) + MOD) % MOD;
	for (int i = k + 1; i <= cnt && pri[i] * pri[i] <= n; ++ i)
		for (int64_t p = pri[i]; p <= n; p *= pri[i])
			ans = (ans + Fp (p) * (S (n / p, i) + (p > pri[i])) % MOD) % MOD;
	return ans;
}

int main () {
	
	std::cin >> N, sqrtn = sqrt (N), Sieve ();
	
	for (int64_t l = 1, r; l <= N; l = r + 1) {
		r = N / (N / l), pos[++ tot] = N / l;
		g1[tot] = S1 (pos[tot]) - 1, g2[tot] = S2 (pos[tot]) - 1;
		pos[tot] <= sqrtn ? id1[pos[tot]] = tot : id2[N / pos[tot]] = tot;
	}
	
	for (int i = 1; i <= cnt; ++ i) {
		for (int j = 1; j <= tot && pri[i] * pri[i] <= pos[j]; ++ j) {
			g1[j] = (g1[j] - pri[i] * (g1[Id (pos[j] / pri[i])] - h1[i - 1]) % MOD + MOD) % MOD;
			g2[j] = (g2[j] - pri[i] * pri[i] % MOD * (g2[Id (pos[j] / pri[i])] - h2[i - 1]) % MOD + MOD) % MOD;
		}
	}
	
	pri[cnt + 1] = 1e18;
	
	std::cout << (S (N, 0) + 1) % MOD << '\n';
	
	return 0;
}