1. 程式人生 > 實用技巧 >Min_25篩 學習小記

Min_25篩 學習小記

前言

為什麼叫學習小記呢?因為暫時除了模板題就沒有做其他的東西了。(霧

這個東西折磨了我一整天,一堆有問題的題解看得我身不如死,只好結合程式碼理解題解,差點死在機房。(話說半天綜合半天競賽真是害人不淺)

為了以後忘了再受荼毒,這裡還是寫一下,如果有人會看到的話,希望可以幫助到吧。(話說這個東西我已經拖了好久了啊!!!)

(話說我怎麼這麼多話說啊?!!)

Min_25 篩

這個東西是由聚聚\(\texttt{Min-25}\)發明了,所以我們稱之為\(\texttt{Min-25}\)篩。(感覺有點民科了)那就不廢話了,進入正言了。

我們考慮對於一個積性函式求和,我們可能會想到杜教篩或者洲閣篩(儘管我不是很會),但是兩種篩法並不是很具有普遍性,特別是杜教篩。

我們考慮一個具體問題,我們現在設一個積性函式\(f(p^k)=p^k(p^k-1)\),其中\(p\)為質數。我們需要求出\(\sum_{i=1}^{n} f(i)\)

我們發現對於這個問題我們可以分成兩個部分進行考慮,一部分為質數,另一部分為合數,所以即為:

\[\sum_{i=1}^{n} f(i)=\sum_{p\in\mathbb{P}\wedge 2\le p\le n} f(p)+\sum_{p\in \mathbb{P}\wedge p^e\le n\wedge p\le \sqrt n} f(p^e)(\sum_{1\le i\le n\wedge \texttt{LPF(i)}> p} f(i)+[e>1]) \]

其中\(\texttt{LPF(i)}\)表示\(i\)的最小質因子(\(\texttt{Lowest Prime Factor}\))。

這個式子很好理解,可能稍微難一點的就是為什麼\(p\le \sqrt n\),這是因為\(\texttt{LPF(n)}\le \sqrt n,s.t. n\notin \mathbb{P}\)

我們發現這個式子似乎不好繼續往下面搞了。

於是,我們構造一種函式\(g(n,i)\)(鬼知道這是怎麼想到的),定義為:

\[g(n,i)=\sum_{x=1}^{n} [x\in \mathbb{P}\vee \texttt{LPF(x)}>\mathbb{P_i}]x^k \]

\(k\)就是我們把\(f(x)\)拆成單項式唯一存在的指數。比如這道題我們是\(f(p)=p(p-1)\),我們就拆成\(p^2-p\),於是\(k\)分別為\(1\)\(2\)

我們考慮如何求出\(g(n,i)\),我們不難得出轉移式:

\[g(n,i)=\left\{\begin{array}{l} g(n,i-1) (\mathbb{P_i}^2>n)\\ g(n,i-1)-\mathbb{P_i}^k(g(\lfloor\frac{n}{\mathbb{P_i}}\rfloor,i-1)-g(\mathbb{P_{i-1}},i-1)) \end{array}\right.\]

這裡還是稍微解釋一下吧。我們的\(g(n,i)\)一定是\(g(n,i-1)\)減去那些是合數且最小質因數大於\(\mathbb{P_{i-1}}\)小於等於\(\mathbb{P_i}\)\(x^k\),因為\(x^k\)是一個完全積性函式,所以我們可以直接拆出\(\mathbb{P_i}^k\)出來。然後自己意會一下應該就可以理解了吧。。。

需要提醒的是,後文裡面的\(g\)其實是把單項式又重新組合成多項式的,如果不是很理解可以見程式碼。

我們發現,我們如果設:

\[S(n,i)=\sum_{x=1}^{n} [\texttt{LPF(x)}>\mathbb{P_i}]f(x) \]

我們就可以在\(S()\)\(g()\)之間連線起某種關係,即:

\[S(n,i)=g(n,|\mathbb{P}|)-\sum_{x=1}^{i-1} f(\mathbb{P_x})+\sum_{\mathbb{P_k}^x\le n\wedge k>i}f(\mathbb{P_k}^x)(S(\lfloor\frac{n}{\mathbb{P_k}}\rfloor,k)+[x>1]) \]

\(|\mathbb{P}|\)就是小於等於\(n\)的質數的個數,還有就是\(\texttt{LPF}\)可以等於自身。這個式子其實就是上面那個式子改了一下,這裡就不贅述了。

我們發現最後的答案就是\(S(n,0)\)

這裡還有一些細節需要提,我們發現這裡面函式所有的自變數都可以表示成\(\lfloor\frac{n}{x}\rfloor\),於是我們其實需要的數只有\(\sqrt n\)個。我們直接離散化一下就好了,空間其實是可以用根號分治做到\(\Theta(\sqrt n)\)的。具體見程式碼。

模板題程式碼

#include <bits/stdc++.h>
using namespace std;

#define Int register int
#define mod 1000000007 
#define ll long long
#define MAXN 200005

int mul (int x,int y){return 1ll * x * y % mod;}
int dec (int x,int y){return x >= y ? x - y : x + mod - y;}
int add (int x,int y){return x + y >= mod ? x + y - mod : x + y;}

ll n,sw,w[MAXN];
int sq,cnt,g1[MAXN],g2[MAXN],id1[MAXN],id2[MAXN],gp1[MAXN],gp2[MAXN],pri[MAXN],vis[MAXN];

int S (ll x,int y){
	if (pri[y] >= x) return 0;
	int p = x <= sq ? id1[x] : id2[n / x],ret = dec (dec (g2[p],g1[p]),dec (gp2[y],gp1[y]));
	for (Int i = y + 1;i <= cnt && 1ll * pri[i] * pri[i] <= x;++ i){
		ll pk = pri[i];
		for (Int k = 1;pk <= x;++ k,pk *= pri[i]){
			int o = pk % mod;ret = add (ret,mul (o,mul (o - 1,S (x / pk,i) + (k != 1))));
		}
	}
	return ret;
}

template <typename T> inline void read (T &t){t = 0;char c = getchar();int f = 1;while (c < '0' || c > '9'){if (c == '-') f = -f;c = getchar();}while (c >= '0' && c <= '9'){t = (t << 3) + (t << 1) + c - '0';c = getchar();} t *= f;}
template <typename T,typename ... Args> inline void read (T &t,Args&... args){read (t);read (args...);}
template <typename T> inline void write (T x){if (x < 0){x = -x;putchar ('-');}if (x > 9) write (x / 10);putchar (x % 10 + '0');}

signed main(){
	read (n);sq = sqrt (n);
	for (Int i = 2;i <= sq;++ i){
		if (!vis[i]) pri[++ cnt] = i,gp1[cnt] = add (gp1[cnt - 1],i),gp2[cnt] = add (gp2[cnt - 1],mul (i,i));
		for (Int j = 1;j <= cnt && i * pri[j] <= sq;++ j){
			vis[i * pri[j]] = 1;
			if (i % pri[j] == 0) continue;
		}
	}
	for (ll l = 1,r;l <= n;l = r + 1){
		r = n / (n / l);w[++ sw] = n / r;g1[sw] = w[sw] % mod;
		g2[sw] = dec (1ll * g1[sw] * (g1[sw] + 1) % mod * (2 * g1[sw] + 1) % mod * 166666668 % mod,1);
		g1[sw] = dec (1ll * g1[sw] * (g1[sw] + 1) % mod * 500000004 % mod,1);
		if (n / r <= sq) id1[n / r] = sw;else id2[r] = sw; 
	}
	for (Int i = 1;i <= cnt;++ i){
		ll sqr = 1ll * pri[i] * pri[i];
		for (Int j = 1;j <= sw && w[j] >= sqr;++ j){
			ll p = w[j] / pri[i];p = (p <= sq ? id1[p] : id2[n / p]);
			g1[j] = dec (g1[j],mul (pri[i],dec (g1[p],gp1[i - 1])));
			g2[j] = dec (g2[j],mul (mul (pri[i],pri[i]),dec (g2[p],gp2[i - 1])));
		}
	}
	write (add (S (n,0),1)),putchar ('\n');
	return 0;
}