「學習筆記」Min25篩

「學習筆記」Min25篩

前言

周指導今天模擬賽五分鐘秒第一題，十分鐘說第二題是 $\text{Min25}​$ 篩板子題，要不是第三題出題人數據範圍給錯了，周指導十五分鐘就 $\text{AK}​$ 了，爲了向 $\text{AK}​$王 學習，真誠的膜拜他，接受紅太陽的指導，下午就學習了一下 $\text{Min25}​$ 篩。

簡介

若是 $f(n)$ 是一個積性函數，且 $f(n)$ 是一個關於 $n$ 的簡單多項式，並能夠快速算出 $f(p^k),\ p\text{ is prime, } k \geq 0$ 的值，那麼大概能夠用 $\text{Min25}$ 篩在 $\mathcal O(\frac{n^\frac{3}{4}}{log_n})$ 求出 $\sum_{i=1}^nf(i)$ 的值。

算法

首先令 $s$ 爲全部 $i\in[1,n], \lfloor\frac{n}{i}\rfloor$ 的集合，並試圖預處理出 $\forall x \in s \sum_{i=2}^x [i\text{ is prime}] f(i)$

可能當 $x$ 不是質數的時候，$f(x)$ 不太好求，咱們先僞裝全部數都是質數，而後用相似埃式篩法的過程把不是質數的 $f(x)$ 值給篩掉。

具體的算法以前，先定義一些東西：

令 $P={prime_1\dots prime_m}$ 表示前 $m$ 個質數的集合，而且 $prime_{m+1}>\sqrt n$ 。

令 $g(n,k) ={\sum_{i=2}^n[i \text{ is prime}\text{ or } minp(i)>prime_k]}f(i) $ 表示前 $n$ 個數中知足是質數或者最小質因子大於第 $k$ 個質數的數僞裝它爲質數求出來的 $f(x)$ 之和。

$g(n, k)​$ 的本質是篩掉了前 $n​$ 個數篩去前 $k​$ 個質數的倍數後剩下的那些數的 $f(x)​$ 之和，顯然，$g(S,m)​$ 就是要預處理的東西。

能夠經過這個東西求 $g$ 了 $$ g(n,k) = g(n,k-1)\ \ \ \ \ \ \text{if } prime_k > \sqrt n \ g(n,k) = g(n,k-1)-f(prime_k)\times [g(\lfloor\frac{n}{prime_k}\rfloor,k-1)-\sum_{i=1}^{k-1} f(prime_i)]\ \ \ \ \ \ \text{if } prime_k \leq \sqrt n $$ 第一個轉移顯然，$prime_k​$ 篩不掉任何數了，第二個轉移考慮把全部在前 $k-1​$ 輪最小質因子不爲 $prime_k​$ 的數貢獻減掉，由於 $f(x)​$ 是積性的直接搞就行了，由於 $g(n,k)​$ 仍然要保留前 $k-1​$ 個質數的貢獻，因此還要加回來。

如今已經預處理出了 $\forall x \in s \sum_{i=2}^x [i\text{ is prime}] f(i)​$ 的值，也就是咱們求出了質數的答案，接下來經過從小到大加入質因子來求出合數的答案。

令 $S(n,k)=\sum_{i=2}^n [minp(i)\geq k] f(i)$ ，這裏的 $f(i)$ 就是真實的值了，不僞裝他是質數，那麼答案就是 $f(1)+S(n,1)$ 。

質數的貢獻以前已經算出來了，就是 $g(n,m)-\sum_{i=1}^{k-1}f(prime_i)$。

考慮枚舉每一個數的最小質因子以及它們冪次獲得合數的轉移 $$ S(n,k)=g(n,m)-\sum_{i=1}^{k-1}f(prime_i)+\sum_{j=k}^m\sum_{t=1}^{prime_{j}^{t+1}\leq \ n}S(\lfloor\frac{n}{prime_j^t}\rfloor,j+1)\times f(prime_{j}^t)+f(prime_j^{t+1}) $$ 其中 $S(\lfloor\frac{n}{prime_j^t}\rfloor,j+1)\times f(prime_{j}^t)$ 算的是後面還有別的質因子的狀況，$f(prime_j^{t+1})$ 算的是後面沒有別的質因子的狀況，能夠感性理解一下。

例題

「Loj6235」 區間素數個數

求 $[1,n]$ 的素數個數，$n \leq10^{11}$

令 $f(x)$ 當 $x$ 是質數的時候爲 $1$ ，$x$ 是其它數的時候爲 $0$ ，雖然不知足這個東西是一個關於 $x$ 的簡單多項式，可是隻求 $f$ 是質數的時候的答案顯然是對的，$g(n,m)$ 就是答案。

code

/*program by mangoyang*/
#include <bits/stdc++.h>
#define inf (0x7f7f7f7f)
#define Max(a, b) ((a) > (b) ? (a) : (b))
#define Min(a, b) ((a) < (b) ? (a) : (b))
typedef long long ll;
using namespace std;
template <class T>
inline void read(T &x){
	int ch = 0, f = 0; x = 0;
	for(; !isdigit(ch); ch = getchar()) if(ch == '-') f = 1;
	for(; isdigit(ch); ch = getchar()) x = x * 10 + ch - 48;
	if(f) x = -x;
}
#define ID(x) ((x) <= T ? id[0][(x)] : id[1][n/(x)])
const int N = 1000005;
ll prime[N], id[2][N], a[N], g[N], n, m, tot, T;
int b[N];
int main(){
	read(n), T = sqrt(n);
	for(ll l = 1; l <= n; l = n / (n / l) + 1){
		a[++m] = n / l, g[m] = a[m] - 1;
		id[a[m]<=T?0:1][a[m]<=T?a[m]:n/a[m]] = m;
	}
	for(int i = 2; i <= T; i++){
		if(!b[i]) prime[++tot] = i;
		for(int j = 1; j <= tot && i * prime[j] <= T; j++){
			b[i*prime[j]] = 1;
			if(i % prime[j] == 0) break;
		}
	}
	for(int i = 1; i <= tot; i++)	
		for(int j = 1; j <= m && prime[i] * prime[i] <= a[j]; j++)
			g[j] -= g[ID(a[j]/prime[i])] - (i - 1);
	cout << g[ID(n)] << endl;
	return 0;
}

<br>

「51NOD1222」 最小公倍數計數

令 $f(n)=\sum_{i=1}^n\sum_{j=i}^n [\text{lcm}(i,j)=n]$，求 $\sum_{i=a}^b f(i),a\leq b\leq 10^{11}$ 。

首先先令 $f'(n)=\sum_{i=1}^n\sum_{j=1}^n[\text{lcm}(i,j)=n]$ ，顯然 $f(n)=\frac{f'(n)+n}{2}$ 。

$$ \sum_{i=1}^n\sum_{j=1}^n[\text{lcm}(i,j)=n]\ =\sum_{i|n}^n\sum_{j|n}^n[\frac{ij}{\gcd(i,j)}=n] \ =\sum_{i|n}^n\sum_{j|n}^n[ij=\gcd(ni,nj)] \ =\sum_{i|n}^n\sum_{j|n}^n[\gcd(\frac{n}{i},\frac{n}{j})=1] \ =\sum_{i|n}^n\sum_{j|n}^n[\gcd(i,j)=1] \ =\sum_{d|n}2^\omega(d) $$

最後一步考慮組合意義，每一種質因子，要麼所有歸 $i$ 要麼所有歸 $j$ ，而後你會發現，$f'(n)$ 是一個積性函數，且知足 $f(prime_i^k)=2k+1$ ，直接 $\text{Min25}$ 篩便可。

code

#include<bits/stdc++.h>
typedef long long ll;
using namespace std;
const int N = 1000005;
namespace Min25{
	#define id(x) ((x) <= T ? id1[x] : id2[n/(x)])
	ll prime[N];
	int b[N], tot, id1[N], id2[N], m;
	ll a[N], g[N], T, n;
	inline void init(){
		m = tot = 0, T = sqrt(n + 0.5);
		for(int i = 2; i <= T; i++){
			if(!b[i]) prime[++tot] = i;
			for(int j = 1; j <= tot && i * prime[j] <= T; j++){
				b[i*prime[j]] = 1;
				if(i % prime[j] == 0) break;
			}
		}
		for(ll l = 1; l <= n; l = n / (n / l) + 1){
			a[++m] = n / l;
			if(a[m] <= T) id1[a[m]] = m; else id2[n/a[m]] = m;
			g[m] = 3ll * (a[m] - 1);
		}
		for(int j = 1; j <= tot; j++)
			for(int i = 1; i <= m && prime[j] * prime[j] <= a[i]; i++){
					g[i] -= g[id(a[i]/prime[j])] - 3ll * (j - 1);
				}
	}
	inline ll solve(ll a, ll b){
		if(a < prime[b]) return 0;
		ll ans = g[id(a)] - 3 * (b - 1);
		for(int j = b; j <= tot && prime[j] * prime[j] <= a; j++)
			for(ll p = prime[j], t = 1; p * prime[j] <= a; t++, p = p * prime[j])
				ans += solve(a / p, j + 1) * (2 * t + 1) + 2 * t + 3;
		return ans;
	}
	inline ll gao(ll x){
		if(x == 0) return 0;
		if(x == 1) return 1;
		n = x, init();
		return (solve(n, 1) + 1 + n) / 2ll;
	}
}
int main(){
	ll a, b;
	cin >> a >> b;
	cout << Min25::gao(b) - Min25::gao(a - 1) << endl;
	return 0;
}