前言

由於其由 Min_25 發明並最早開始使用，故稱「Min_25 篩」。

• 從此種篩法的思想方法來說，其又被稱為「Extended Eratosthenes Sieve」。

其可以在 \(O(\frac{n^{\frac{3}{4}}}{logn})\) 的時間複雜度下解決一類 積性函式 的字首和問題。
要求： \(f(p)\) 是一個關於 \(p\) 的項數較少的多項式或可以快速求值； \(f(p^{c})\) 可以快速求值。
摘自 OI wiki

解析

例題

積性函式 \(f(p^k) = p^k(p^k - 1)\)
求：\(\sum_{i = 1} ^ n f(i)\)

準備 \(1\)

規定 \(p\) 代表質數, \(minp(i)\) 代表 \(i\) 的最小質因子（其中 \(minp(p) = p\)），則：
\(\sum_{i = 1} ^ n f(i) = \sum_{p \leq n} f(p) + \sum_{i \leq n \ and \ i \notin p} f(i)\)
這一步表示將答案拆成質數的答案和合數的答案。
那麼對於合數部分：
根據積性函式的性質，則有：
\(f(i) = \prod_{k = 1} ^ m f(p_k ^ {c_k})\)。
提盡 \(i\) 的最小質因子 \(p_e\)，將 \(i\) 分成兩部分：
\(p_e ^ {c_e} \times \frac{i}{p_e ^ {c_e}}\)

• \(p\) 是否是 \(i\) 的最小質因子
• \(p\) 是否提盡

這兩個條件均可通過 \(minp(\frac{i}{p ^ e})\) 是否大於 \(p\) 來判斷：
對於 \(minp(\frac{i}{p ^ e}) < p\) 的 \(i\)，那麼不滿足條件 \(1\)。
對於 \(minp(\frac{i}{p ^ e}) = p\) 的 \(i\)，那麼不滿足條件 \(2\)。
這些均可通過對 \(i\) 分解質因數來理解。
那麼對於 \(f(p ^ e)\) 由題意可以直接求，對於 \(f(\frac{i}{p ^ e})\) 則是同樣的求法，直到 \(\frac{i}{p ^ e}\) 是 \(p_k ^ {c_k}\) 可以直接求。
而上述的約束則保證了不會多求漏求。

準備 \(2\)

所以我們現在需要求出所有 \(f(p ^ k)\) 的值。
因為 \(n\) 的範圍很大，所以我們沒法線性篩求出 \(f\) 的字首和。
考慮 DP。。。
設 \(h_k(i) = i ^ k\)，顯然 \(h_k\) 是一個完全積性函式。
設 \(g_k(n,j) = \sum_{i = 1} ^ n [i \in Prime \ or \ minp(i) > p_j] h_k(i)\)，其中 \(p_j\) 表示第 \(j\) 個質數。
換句話說 \(g_k(n,j)\) 也就是 \(1\) ~ \(n\) 中質數或者最小質因數大於 \(p_j\) 的 \(i\) 的 \(h_k(i)\) 的字首和。
\(g_k\) 的作用在下文，但也可以先拋開問題本身，看懂 \(g_k\) 的求法。
考慮轉移。。。
對於 \(g_k(n, j - 1)\) 根據定義，其中包含兩部分：

• \(g_k(n, j)\)
• 最小質因數等於 \(p_j\) 的 \(h_k(i)\) 的和。

所以可由 \(g_k(n, j - 1)\) 減去多餘的部分得到 \(g_k(n, j)\)。
所以有狀態轉移方程如下：
\(g_k(n,j) = g_k(n, j - 1) - p_{j} ^ k (g(\frac{n}{p_{j} ^ k}, j - 1) - g(p_j, j - 1))\)
從狀態轉移方程可以看出：
最小質因數等於 \(p_j\) 的 \(h_k(i)\) 的和即為 \(g_k(n, j - 1) - p_{j} ^ k (g(\frac{n}{p_j}, j - 1) - g(p_{j - 1}, j - 1))\)。
原因如下：
根據定義不難得到 \(g_k(p_j, j - 1)\) 即為前 \(j - 1\) 個質數的 \(k\) 次方和。
而對於 \(g_k(\frac{n}{p_j}, j - 1)\)，包含了兩部分：

• \(1\) ~ \(p_{j - 1}\) 的質數和
• 最小質因數大於等於 \(p_j\) 的合數和，及 \(p_{j}\) ~ \(\frac{n}{p_j}\) 的質數和

注意第二點，如果我們將這些最小質因數大於等於 \(p_j\) 的合數或者後一部分的質數同時乘上 \(p_j\) 那麼我們得到的是不是都是最小質因數等於 \(p_j\) 的合數，這點應該很好理解。
那也就是說得到了我們想要篩去的合數。
因為我們只知道上述兩部分的和 \(g_k(\frac{n}{p_j}, j - 1)\)，所以在乘的過程中第一部分也會被乘上一個 \(p_j\)，那麼我們應該把這一部分質數減去。
也就是 \(g_k(p_j, j - 1)\)。
又因為 \(h_k\) 是完全積性函式，而我們談論的和是 \(h_k\) 的累加，根據乘法分配律，我們可以對 \(g_k\) 運用完全積性函式的性質（\(f(ab) = f(a)f(b)\)）。
那麼要得到最小質因數等於 \(p_j\) 的 \(h_k(i)\) 的和，我們可以直接給我們的第二部分乘上 \(h_k(p_j)\) 的值，也就是 \(p_j ^ k\)，得到我們要求的函式值。
所以關於狀態轉移方程的解釋就結束了。
關於初始化：
顯然對於 \(g_k(n, 0) = \sum_{i = 1} i ^ k\)

求解

而我們得到的 \(g_k\) 有什麼用呢？
設 \(g(n, j) = \sum_{i = 1} ^ n [i \in Prime \ or \ minp(i) > p_j] f(i)\)。
因為 \(f(p) = p(p - 1) = p ^ 2 - p\)。
那麼將滿足條件的 \(f(i)\) 累加，就得到兩部分：

• \(i ^ 2\) 的和
• \(i\) 的和的相反數

如果令上文的 \(k\) 分別等於 \(1, 2\)。
那麼就有： \(g(n, j) = g_2(n, j) - g_1(n, j)\)。
如果設 \(p_x\) 為小於等於 \(\sqrt{n}\) 的最大的質數，那麼 \(1\) ~ \(n\) 的質數的 \(f\) 的和 \(g(n, x)\) 記為 \(g(n)\)。
設 \(S(n, j) = \sum_{minp(i) > p_j} f(i)\)。
顯然我們也需要 DP 。。。
同樣的，考慮 \(S(n, j)\) 的組成：

• 質數部分：大於 \(p_j\) 小於等於 \(n\) 的質數
• 合數部分

設 \(sum_k(x)\) 表示小於等於 \(p_x\) 的質數的 \(k\) 次方的和。
那麼質數部分就等於：
\(g(n) - （sum_2(j) - sum_1(j)）\)
對於合數部分：
我們模仿 \(g_k\) 的求法
因為 \(f\) 不是完全積性函式，而是積性函式，我們不能像求 \(g_k\) 那樣只提一個 \(p_j\)。
所以我們要把 \(p_j\) 提盡。
得到：
\(\sum_{k > j \ and \ p_k ^ e \leq n} f(p_k ^ e)(S(\frac{n}{p_k ^ e} + [e \neq 1]))\)
利用 \(S(\frac{n}{p_k ^ e}\) 的思想和求 \(g_k\) 的時候完全一樣，這裡就不在贅述了，實在不懂的同學可以手推一下。
後面加上一個 \([e \neq 1]\) 是因為如果 \(e > 1\)，那麼\(p ^ e\) 在我們提 \(p_k ^ e\) 時沒有算進去，而 \(e = 1\) 時，在質數部分就算了，就不用了算進去。

下標的離散化

因為 \(n = 1e10\)，直接開肯定不行。
由性質 \(\lfloor \frac{\lfloor \frac{x}{a} \rfloor}{b} \rfloor = \lfloor \frac{x}{ab} \rfloor\)，在遞迴的時候，無論我們把 \(n\) 除以幾，得到的都是 \(n\) 除以某個數，所以我們可以直接只儲存 \(n / x\) 的值。
來自 wucstdio 大佬

程式碼

#include<cstdio>
#include<cmath>
#include<cstring>

using namespace std;

typedef long long LL;

const LL mod = 1e9 + 7, inv2 = 500000004, inv3 = 333333336;
const int N = 1e5 + 5;

LL n, Prime[N + 5], num, sp1[N + 5], sp2[N + 5], tot =  0, w[3 * N], g1[3 * N], g2[3 * N], sqr;
int id1[N + 5], id2[N + 5];
bool notPrime[N + 5];

void init() {
	memset(notPrime, false, sizeof(notPrime));
	num = 0;
	for(int i = 2; i <= N; i++) {
		if(!notPrime[i]) {
			Prime[++num] = i;
			sp1[num] = (sp1[num - 1] + i) % mod;
			sp2[num] = (sp2[num - 1] + 1ll * i * i % mod) % mod;
		}
		for(int j = 1; j <= num && Prime[j] * i <= N; j++) {
			notPrime[i * Prime[j]] = true;
			if(i % Prime[j] == 0) break;
		}
	}
}

LL s(LL x, int y) {
	if(Prime[y] >= x) return 0;
	LL k = x <= sqr ? id1[x] : id2[n / x];
	LL ans = (g2[k] - g1[k] + mod - (sp2[y] - sp1[y]) + mod) % mod;
	for(int i = y + 1; i <= num && Prime[i] * Prime[i] <= x; i++) {
		LL P = Prime[i];
		for(int e = 1; P <= x; e++, P = P * Prime[i]) {
			LL xx = P % mod;
			ans = (ans + xx * (xx - 1) % mod * (s(x / P, i) + (e != 1))) % mod;
		}
	}
	return ans % mod;
}

int main() {
	init();
	scanf("%lld", &n);
	sqr = sqrt(n);
	for(LL l = 1, r; l <= n; l = r + 1) {
		r = (n / (n / l));
		w[++tot] = n / l;
		g1[tot] = w[tot] % mod;
		g2[tot] = g1[tot] * (g1[tot] + 1) / 2 % mod * (2 * g1[tot] + 1) % mod * inv3 % mod - 1;
		if(g2[tot] < 0) g2[tot] += mod;
		g1[tot] = g1[tot] * (g1[tot] + 1) / 2 % mod - 1;
		if(g1[tot] < 0) g1[tot] += mod;
		if(w[tot] <= sqr) id1[w[tot]] = tot;
			else id2[n / w[tot]] = tot;
	}
	for(int i = 1; i <= num; i++)
		for(int j = 1; j <= tot && Prime[i] * Prime[i] <= w[j]; j++) {
			LL k = w[j] / Prime[i];
			k = k <= sqr ? id1[k] : id2[n / k];
			g1[j] -= Prime[i] * (g1[k] - sp1[i - 1] + mod) % mod;
			g2[j] -= Prime[i] * Prime[i] % mod * (g2[k] - sp2[i - 1] + mod) % mod;
			g1[j] %= mod, g2[j] %= mod;
			if(g1[j] < 0) g1[j] += mod;
			if(g2[j] < 0) g2[j] += mod;
		}
	printf("%lld", (s(n, 0) + 1) % mod); 
	return 0;
}