Min_25 篩學習筆記

\(\text{Min}\_25\) 篩學習筆記

事實上我又學習了一個有點春的篩法。\(\text{Min}\_25\) 篩用於求解積性函式的字首和，即形如 \(g(n)=\sum_{i=1}^{n}f(i)\) 形式的函式 \(g\)。

眾所周知，樸素篩法之所以無法做到低於線性是因為列舉了區間內的每一個數，那麼我們想要做到低於線性，就必然需要做到透過一種方法將所有數字分成兩類，用過一類的求解輔助另一類的求解，不難想到質數和合數。那麼我們分開考慮如何求解質數和合數。對於 \(1\) 這個既不是質數也不是合數的數，我們選擇將其的貢獻最後單獨統計，畢竟一個積性函式要麼所有位置都是 \(0\)，要麼 \(f(1)=1\)。

質數求解

首先考慮如何求解質數的貢獻。眾所周知，冪函式是一種完全積性函式，假若當 \(x=p\) 時 \(f(x)=\sum x^i\)，那麼我們便可以將每一項分開求解。更具體的我們嘗試統計 \(\sum_{i\in P}i^k\) 的值。

埃篩是一種複雜度略差於線性篩的篩法，但只表現為擁有較大的常數，在這裡我們借鑑埃篩的思路。考慮埃篩利用每一個範圍內的質數進行過濾，更具體的，對於每一個質數 \(p\)，其都將 \(\ge p^2\) 且\(\bmod p=0\) 的數篩掉。那麼我們考慮所有 \(\le n\) 的合數一定有一個 \(\le\sqrt{n}\) 的質因子，因此我們只需要篩出前 \(\sqrt{n}\) 個素數即可。但是埃篩的時間複雜度依舊不可接受。考慮用 \(\text{dp}\) 最佳化這個過程，我們設 \(g_{i,j}\) 表示 \([1,i]\) 之間的數在經過 \(j\) 輪篩之後所剩餘的數的貢獻，那麼我們需要知道的就是 \(g_{n,|P_{\sqrt n}|}\)。首先有初始值

\[g_{n,0}=\sum_{i=2}^{n}i^k \]

那麼考慮因為質數 \(p_j\) 只會篩掉 \(\ge p_j^2\) 的數，那麼我們可以推斷出轉移方程

\[g_{i,j}=g_{i,j-1}(i<p_j^2) \]

接著考慮 \(i\ge p_j^2\) 的時候，因為完全積性函式的性質，所以我們單獨將被篩掉的合數的質因子 \(p_j\) 提出來。那麼剩下的部分應該是最小質因子 \(\ge p_j\) 的部分，因為如果最小質因子 \(<p_j\) 一定會被之前的某個素數篩掉，我們可以用 \(g_{\lfloor\frac{i}{p_j},j-1\rfloor}\) 來表示，但是這裡麵包含了 \([1,p_{j-1}]\) 的質數，不應該減去這部分，因此還需要將這部分貢獻加回來，所以有轉移方程

\[g_{i,j}=g_{i,j-1}-p_j^k(g_{\lfloor\frac{i}{p_j},j-1\rfloor}-g_{p_{j-1},j-1})(i\ge p_j^2) \]

不難發現可以用滾動陣列最佳化掉第二維，但是時空複雜度依舊沒有消下去。首先後半部分其實沒有必要利用 \(g\) 陣列來表示，反倒可以利用素數字首和陣列 \(s\)，而這個陣列後面也會遇到。接著考慮可能被轉移的位置 \(\lfloor\frac{i}{p_j}\rfloor\)，如果從搜尋的角度看，就是 \(n,\lfloor\frac{n}{p_1}\rfloor,\lfloor\frac{\lfloor\frac{n}{p_1}\rfloor}{p_2}\rfloor,\cdots\)，有一個化簡技巧就是 \(\lfloor\frac{\lfloor\frac{n}{a}\rfloor}{b}\rfloor=\lfloor\frac{n}{ab}\rfloor\)。那麼我們發現有用的下標其實只有 \(\lfloor\frac{n}{a}\rfloor\)，那麼我們可以利用數論分塊在 \(O(\sqrt{n})\) 的時間複雜度內處理出所有可能的值。

合數求解

在求解完素數之後我們就要醞釀著統計合數對答案的貢獻了。因為給定的函式是積性函式的緣故，因此我們可以將合數表示為最小質因子的次方和另一個數相乘的結果，因為後半部分的最小質因子不能小於原來的最小質因子，因此我們令 \(S(n,b)\) 表示範圍 \([1,n]\) 內所有最小質因子 \(>p_b\) 的數的貢獻和，那麼所求即為 \(S(n,0)\)。遞迴求解這個東西。當 \(n\le p_b\) 時，直接返回 \(0\) 即可。否則我們先統計所有素數的答案，然後進行最小質因子和其次冪的列舉，並將對應的答案貢獻累加，更具體的，有轉移方程

\[S(n,b)=g_{n,|P_{\sqrt{n}}|}-\sum_{i=1}^{b}f(p_i)+\sum_{i=b+1}^{|P_{\sqrt{n}}|}\sum_{1\le j,p_i^{j}\le n}(S(\lfloor\frac{n}{p_i^{j}}\rfloor,i)+[j\ne1]f(p_i^j)) \]

我們發現中間的求和可以用之前的質數字首和快速求解，而剩餘部分自行遞迴即可。整體的複雜度是 \(O(\frac{n^{\frac{3}{4}}}{\log n})\)。

\(\text{Min}\_25\) 篩的使用場景包括但不限於所有可以快速求出質數及其冪次處的函式值，並且其在質數處的取值可以視作關於 \(p\) 的多項式。

事實上在程式碼實現時，有一些不加不可的最佳化，和一些美化程式碼的最佳化，只能說自行選擇了。程式碼以 \(f(x)=x\) 作為演示

程式碼

#include <bits/stdc++.h>
using namespace std;
#define ll long long
#define Get(x) (x<=S?id1[x]:id2[n/(x)])
#define f(p) (p)

const int N=1e5+5,P=1e9+7;

int S,cnt,tot;
int p[N],id1[N],id2[N];
ll n;
ll sum[N],g[N<<1],val[N<<1];//n/i 的值最多有 2sqrt(n) 個，結論不要記錯了
bool isp[N];

void Euler_S(){
	isp[1]=true;
	for(int i=2;i<=S;i++){
		if(!isp[i]){
			p[++cnt]=i;
			sum[cnt]=(sum[cnt-1]+i)%P;
		}
		for(int j=1;p[j]<=S/i;j++){
			isp[i*p[j]]=true;
			if(i%p[j]==0)continue;
		}
	}
	return ;
}

ll SumF(ll a,int b){
	if(a<=p[b])return 0;
	int id=Get(a);
    ll res=g[id]-sum[b];
	for(int i=b+1;i<=cnt&&p[i]<=a/p[i];i++){//注意這裡後半部分的判斷是必加的，不加也能夠得出正確答案，但因為常數過大被卡爆了
		ll pk=p[i];
		for(int j=1;pk<=a/p[i];j++,pk*=p[i])res=(res+f(pk)*SumF(a/pk,i)+f(pk*p[i]))%P;
        //這部分和上面所說的轉移方程略有不同，因為考慮到如果 pk*p>a，那麼呼叫過去答案一定為 0，那麼就沒有必要進行這次操作
        //並且後半部分變成 f(p^(j+1))，這樣既能夠避免重複統計質數，也能剛好覆蓋到最大的數，同時能夠保障 pk 不超過 ll 範圍
        //這部分最佳化可以自行選擇
	}
	return res;
}

int main(){
	ios::sync_with_stdio(0);
	cin.tie(0),cout.tie(0);
	cin>>n;
	S=sqrt(n);
	Euler_S();
	for(ll l=1,r,v;l<=n;l=r+1){
		r=n/(v=n/l);
		if(v<=S)id1[v]=++tot;
		else id2[n/v]=++tot;//兩種分開儲存
		val[tot]=v;
		v%=P;
		g[tot]=v*(v+1)/2%P-1;//別忘了初始值沒有包括 1
	}
	for(int j=1;j<=cnt;j++){
		for(int i=1;i<=tot&&p[j]<=val[i]/p[j];i++){
			int id=Get(val[i]/p[j]);
			g[i]=(g[i]-(g[id]-sum[j-1]))%P;//直接用字首和代替 g 陣列
		}
	}
	cout<<(SumF(n,0)+1+P)%P;//記得單獨加上 1 的貢獻
	return 0;
}