Min_25 篩


Min_25 篩是一種亞線性篩法,可以在 \(\mathcal{O}(\frac{n^{\frac{3}{4}}}{\log n})\) 的時間復雜度下快速算出形如:

\[\sum_{i=1}^n f(i) \]

的值,不過一般比較好實現的方法被證明復雜度是 \(\mathcal{O}(n^{1-\epsilon})\)

本文約定 \(\mathbb{P}\) 表示質數集,\(\mathbb{P}'\) 表示所有 \(\le \sqrt{n}\) 的質數集。

對於函數的限制

在了解一個篩法的具體內容前,先了解其適用范圍是有助於推導過程時的理解和應用時不至於用錯的。對於 Min_25 篩,其要求求和的 \(f(i)\) 滿足如下性質:

  • 積性函數,即對於任意的 \(\gcd(x,y)=1\),都有 \(f(x)f(y)=f(xy)\)
  • \(f(p),p\in\mathbb{P}\) 能被表示為項數比較小的多項式。
  • \(f(p^c),p\in\mathbb{P}\) 能夠快速求值。

接下來,我們就找一個滿足以上條件的函數並以它為例子進行介紹吧。

P5325 【模板】Min_25篩

定義積性函數 \(f(x)\),其中 \(f(p^k)=p^k(p^k-1),p\in\mathbb{P}\),求:

\[\sum_{i=1}^nf(i)\bmod{10^9+7} \]

(\(1\le n\le 10^{10}\))

對於和式的恰當拆分

首先我們注意到許多性質是針對質數來的,所以我們考慮把質數從求和式中單獨拉出來:

\[\sum_{i=1}^nf(i)=\sum_{p=1}^n[p\in\mathbb{P}]f(p)+\sum_{i=1}^nf(i)[i\notin\mathbb{P}\operatorname{or}i=1] \]

但對於后面的和式怎么辦呢?我們想到因為是積性函數,所以合數可以提取出來質因數與質數扯上關系,枚舉合數的 \(\rm LPF\)\(\rm lowest\ prime\ factor\),最小質因數)及其次數有:

\[\sum_{p=1}^n[p\in\mathbb{P}]f(p)+\sum_{p=1}^{\sqrt{n}}[1\le p^e\le n\operatorname{and}p\in\mathbb{P}]f(p^e)\left(\sum_{i=1}^{\lfloor\frac{n}{p^e}\rfloor}[\mathrm{LPF}_i>p]f(i)\right) \]

提取出來 \(p^e\) 后如果原數滿足 \(\mathrm{LPF}=p\) 則剩下的數應該滿足 \(\mathrm{LPF}_i>p\)

引入遞推思想

(不知道 Min_25 神仙咋想到的啊Orz)考慮遞推:設 \(g(n,j)\) 表示:

\[\sum_{i=1}^n i^k[i\in\mathbb{P}\operatorname{or}\mathrm{LPF}_i>p_j] \]

這里 \(i^k\) 是把 \(f(p)\) 作為多項式拆分成單項式的結果,所以如果 \(f(p)\)\(l\) 項,那我們就要求 \(l\)\(g(n,j)\),並用它們的線性組合組合出 \(f(p),p\in\mathbb{P}\)

\(g(n,j)\) 的定義說人話就是對於要么是質數,要么 \(\rm LPF\) 大於 \(p_j\) 的數求 \(k\) 次方和。根據這個,我們能發現 \(\sum_{p=1}^nf(p)[p\in\mathbb{P}]\) 就能用 \(g(n,|\mathbb{P}'|)\) 的線性組合來表示了,其中我們用 \(\mathbb{P}'\) 是因為一個合數 \(d\)\(\rm LPF\) 一定 \(\le \sqrt{d}\)

接下來的問題就是為 \(g(n,j)\) 找到一個遞歸式,也即如何實現 \(g(n,j-1)\rightarrow g(n,j)\) 的轉移。考慮這個轉移的實質其實是把限制拉大了,所以一定是去掉一些數,而去掉的數應該滿足 \(\mathrm{LPF}=p_j\)。既然滿足 \(\mathrm{LPF}=p_j\),我們就試着提取一個 \(p_j\) 出來:

\[g(n,j)=g(n,j-1)-p_j^k\left(g\left(\dfrac{n}{p_j},j-1\right)-g(p_{j-1},j-1)\right) \]

解釋一下這個式子。首先要發現,\(i^k\) 是完全積性函數,即 \(\forall x,y\) 都有 \(f(x)f(y)=f(xy)\),所以即使 \(\dfrac{n}{p_j}\) 不一定與 \(p_j\) 互質也能提取(而上文 \(f\) 僅為積性函數,所以必須一口氣提取完所有的 \(p\) 才能保證互質)。然后分別看一下括號里面的兩項。\(g\left(\left\lfloor\dfrac{n}{p_j}\right\rfloor,j-1\right)\) 關於的是提取完 \(p_j\) 后的數,而一個數如果想 \(\mathrm{LPF}=p_j\),提取完 \(p_j\) 后其 \(\rm LPF\) 應該 \(> p_{j-1}\)。但這有個問題,\(\le p_{j-1}\) 的質數會被統計到這里面,但是它們並不滿足 \(\mathrm{LPF}>p_{j-1}\) 所以要減去。\(g(p_{j-1},j-1)\) 就是干這個的,顯然 \(\le p_{j-1}\) 的數不可能滿足 \(\mathrm{LPF}>p_{j-1}\),所以這個統計的就是 \(\le p_{j-1}\) 的質數。

求遞歸式

我們剛剛求出了關於 \(g(n,j)\) 的遞歸式,但發現求 \(g(n,j)\) 還得用到 \(g\left(\left\lfloor\frac{n}{p_j}\right\rfloor,j-1\right)\),求 \(g\left(\left\lfloor\frac{n}{p_j}\right\rfloor,j-1\right)\) 又得用到 \(g\left(\left\lfloor\frac{\left\lfloor\frac{n}{p_j}\right\rfloor}{p_k}\right\rfloor,j-2\right)\) 等等。以此類推,難道我們對於所有 \(1\sim n\) 都求出 \(g(n,j)\) 嗎?顯然復雜度不能接受。但我們發現上述過程中出現了一個很熟悉的結構,有結論:

\[\left\lfloor\dfrac{\left\lfloor\dfrac{n}{a}\right\rfloor}{b}\right\rfloor=\left\lfloor\dfrac{n}{ab}\right\rfloor \]

那容易發現,不管遞推多少次,我們需要的值總能被表示為 \(g\left(\left\lfloor\dfrac{n}{x}\right\rfloor,j\right),x\in\mathbb{N}_+\) 的形式,我們只需要對 \(\left\lfloor\dfrac{n}{x}\right\rfloor\)\(g\) 就好了,而它僅有 \(\mathcal{O}(\sqrt{n})\) 種,可以接受。

有一個實現上的小問題,我們在套上述遞歸式計算的時候,肯定要給每個 \(\left\lfloor\dfrac{n}{x}\right\rfloor\) 一個編號,但在遞歸的時候我們需要實現從 \(\left\lfloor\dfrac{n}{x}\right\rfloor\) 到其編號的轉化,而前者的值域是 \([0,n]\),如果用 std::map 的話復雜度會多一個 \(\log\)。不過也不難實現 \(\mathcal{O}(1)\) 的轉化,考慮對於所有 \(\le \sqrt{n}\)\(\left\lfloor\dfrac{n}{x}\right\rfloor\) 開一個桶記錄 \(\left\lfloor\dfrac{n}{x}\right\rfloor\rightarrow \mathrm{index}_{\lfloor\frac{n}{x}\rfloor}\) 的關系,對於 \(>\sqrt{n}\)\(\left\lfloor\dfrac{n}{x}\right\rfloor\) 用另一個桶記錄 \(\left\lfloor\dfrac{n}{\left\lfloor\dfrac{n}{x}\right\rfloor}\right\rfloor\rightarrow \mathrm{index}_{\lfloor\frac{n}{x}\rfloor}\) 的關系。這樣鍵值的值域就保持在 \([0,\sqrt{n}]\) 了,可以數組存。

另一個遞歸式

考慮設 \(S(n,x)\) 表示:

\[\sum_{i=1}^n f(i)[\mathrm{LPF}_i>p_x] \]

說人話就是所有 \(\mathrm{LPF}>p_x\) 的數求和,則顯然答案就是 \(S(n,0)+1\)(因為一開始求 \(g\) 的時候把 \(1\) 排除在外了)。接下來考慮怎么計算它。注意到符合條件的數有兩種:

  • \(> p_x\) 的質數。
  • \(\mathrm{LPF}_i>p_x\) 的合數。

首先是 \(> p_x\) 的質數,這可以通過我們剛剛算出的 \(g\) 求,即:

\[g(n,|\mathbb{P}'|)-sp_x \]

其中 \(sp_x\) 表示 \(\sum_{p=1}^{p_x}p^k[p\in\mathbb{P}]\),可以通過線性篩在 \(\mathcal{O}(\sqrt{n})\) 的時間前綴和出來。

然后是 \(\mathrm{LPF}_i>p_x\) 的合數,這時候我們考慮枚舉 \(\rm LPF\) 及其次數,得到:

\[\sum_{p_k^e\le n\operatorname{and}k>x}f(p_k^e)\left(S\left(\dfrac{n}{p^e_k},k\right)+[e\ne 1]\right) \]

這里的提取跟最開始拆分 \(g\) 的思路類似,只不過這次因為不是完全積性函數是全部提取完的,\(\rm LPF\) 應該要 \(>p_k\)。然后是 \([e\ne 1]\),顯然 \(f(p_k^e)\) 應該算上,但這樣提取是不會被 \(S\) 包含的,應該 \(+1\)\(f(p_k^e)\) 包含進去。但當 \(e=1\) 的時候,\(f(p_k^e)\) 是質數,已經被算過了,不能加上。

這樣把上面倆式子加一起就可以算了,這里采用遞歸的形式實現最方便,由於某個我並不知道的定理,可以不用記憶化也能保證復雜度。最后

總結

Min_25 篩的基本步驟:

  1. 將待求函數在質數情況下的表達式轉化成 \(\sum i^k\) 的形式。
  2. 線性篩出 \(1\sim \sqrt{n}\) 之間的質數,並算出其 \(k\) 次方前綴和。
  3. 對於形如 \(\lfloor\frac{n}{x}\rfloor\) 的數算出 \(g(\lfloor\frac{n}{x}\rfloor,0)\),注意去掉 \(1^k\)
  4. 套遞歸式算 \(g(n,j)\),這里 \(g\) 可以用滾動數組.
  5. 遞歸計算 \(S(n,x)\),無需記憶化。答案即為 \(S(n,0)+1\)

板子題代碼:

#include <cmath>
#include <cstdio>
const int N = 1e6 + 10, mod = 1e9 + 7, inv3 = 333333336, inv2 = 500000004; 
typedef long long ll; ll n, w[N]; int p[N], vis[N], sp1[N], sp2[N], tp, sqrn;
// 線性篩出 <= \sqrt{n} 的質數,並求出 p^k 的前綴和
inline void getP(int m)
{
	vis[1] = 1;
	for (int i = 2; i <= m; ++i)
	{
		if (!vis[i]) p[++tp] = i;
		for (int j = 1; j <= tp && (ll)p[j] * i <= m; ++j)
		{
			vis[i * p[j]] = 1;
			if (i % p[j] == 0) break;
		}
	}
	for (int i = 1; i <= tp; ++i) sp1[i] = (sp1[i - 1] + p[i]) % mod;
	for (int i = 1; i <= tp; ++i) sp2[i] = (sp2[i - 1] + (ll)p[i] * p[i] % mod) % mod;
}
int ind1[N], ind2[N], tot, g1[N], g2[N]; 
int S(ll x, int y)
{
	if (p[y] >= x) return 0;
	int k = x <= sqrn ? ind1[x] : ind2[n / x];
	int ans = ((ll)g2[k] - g1[k] - sp2[y] + sp1[y] + mod + mod) % mod;
	for (int i = y + 1; i <= tp && (ll)p[i] * p[i] <= x; ++i)
	{
		ll pe = p[i];
		for (int e = 1; pe <= x; ++e, pe *= p[i])
		{
			int xx = pe % mod;
			ans = (ans + (ll)xx * (xx - 1) % mod * (S(x / pe, i) + (e != 1)) % mod) % mod;
		}
	}
	return ans;
}
int main()
{
	scanf("%lld", &n); sqrn = sqrt(n); getP(sqrn);
	// 求出 g(\lfloor\frac{n}{x}\rfloor,0) 的值,並存一下 \lfloor\frac{n}{x}\rfloor 對應的下標
	for (ll l = 1, r, tn; l <= n; l = r + 1)
	{
		r = n / (n / l); w[++tot] = n / l; tn = w[tot] % mod;
		g1[tot] = tn * (tn + 1) % mod * inv2 % mod - 1;
		g2[tot] = tn * (tn + 1) % mod * (2 * tn + 1) % mod * inv2 % mod * inv3 % mod - 1;
		if (w[tot] <= sqrn) ind1[w[tot]] = tot;
		else ind2[n / w[tot]] = tot;
	}
	// 根據遞歸式算出 g(n,j) 的值
	for (int i = 1; i <= tp; ++i)
		for (int j = 1; j <= tot && (ll)p[i] * p[i] <= w[j]; ++j)
		{
			int k = (w[j] / p[i]) <= sqrn ? ind1[w[j] / p[i]] : ind2[n / (w[j] / p[i])];
			(g1[j] -= (ll)p[i] * (g1[k] - sp1[i - 1] + mod) % mod) %= mod;
			(g2[j] -= (ll)p[i] * p[i] % mod * (g2[k] - sp2[i - 1] + mod) % mod) %= mod;
			(g1[j] += mod) %= mod; (g2[j] += mod) %= mod;
		}
	// 遞歸求 S(n,0),根據某個定理,不需要記憶化也能保證復雜度
	printf("%d\n", (S(n, 0) + 1) % mod); return 0;
}


免責聲明!

本站轉載的文章為個人學習借鑒使用,本站對版權不負任何法律責任。如果侵犯了您的隱私權益,請聯系本站郵箱yoyou2525@163.com刪除。



 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM