min25篩學習筆記

話說我們現在要求一個函數\(f\)的前綴和。即求\(F(n)=\sum_{i=1}^nf(i)\)。

min25篩這個算法的主要思想是把\(1...n\)這些數按質數和合數分類，然后分別考慮質數和合數的貢獻。

STEP1 質數貢獻

我們嘗試先解決一個小問題：求\(G(m)=\sum_{i=1}^m[i\in prime]f(i)\)，即\(m\)以下的質數的\(f\)值之和。其中\(m\in\{\lfloor\frac{n}{1}\rfloor,\lfloor\frac{n}{2}\rfloor,...,\lfloor\frac{n}{n}\rfloor\}\)，共有\(O(\sqrt{n})\)種取值。

這玩意很難直接求，於是考慮DP。設\(g(j,m)=\sum_{i=1}^m[i\in prime \; or \; minp(i)>p[j]]f'(i)\)，其中\(minp(i)\)表示\(i\)的最小質因子。這個式子的含義是“\(i\)為質數或\(i\)的最小質因子\(>p[j]\)的函數值\(f'(i)\)”之和。其中\(f'(i)\)並不一定是我們要求的\(f(i)\)，它可以是我們構造出來的另一個函數，但它必須滿足以下條件：

1. \(f'\)在質數處的取值與\(f\)相同。即\(\forall p\in prime\)有\(f'(p)=f(p)\)。
2. \(f'\)是完全積性函數。
3. \(f'\)可以快速求前綴和。

顯然DP的邊界是\(g(0,m)=\sum_{i=2}^mf'(i)\)。

考慮轉移，求\(g(j,m)\)。當\(p[j]^2>m\)時，\(m\)以內不可能有任何最小質因子是\(p[j]\)的合數，即\(p[j]\)無法篩掉任何數，因此\(g(j,m)=g(j-1,m)\quad (p[j]^2>m)\)。

否則，我們只需要在\(g(j-1,m)\)里“篩掉”\(\leq m\)的最小質因子為\(p[j]\)的合數即可。而它們都可以表示為“\(p[j]\times\)一個\(\lfloor\frac{m}{p[j]}\rfloor\)以內的最小質因子大於等於\(p[j]\)的數”。於是有：

\[g(j,m)=g(j-1,m)-f'(p[j])\times \big(g(j-1,\lfloor\frac{m}{p[j]}\rfloor)-\sum_{i=1}^{j-1}f'(p[i])\big)\qquad (p[j]^2\leq m) \]

這時\(f'\)作為一個完全積性函數的好處就體現出來了，它可以直接乘以后面的一坨東西而不用擔心是否互質。請讀者認真理解這個式子的含義。

令\(|P|\)表示\(\leq m\)的質數個數，則\(g(|P|,m)\)就是我們前面要求的\(G(m)\)的值了。顯然，我們可以把DP的第一維滾掉。

在具體實現時，我們先預處理出\(\sqrt{n}\)以內的質數（線性篩），以及所有\(m\)的值（用數論分塊）。暴力的做法我們可以用STL-map來記錄所有\(m\)的值的編號，但這里可以更巧妙一些。我們開兩個數組id1和id2，大小都是\(\sqrt{n}\)，巧妙地利用“對所有\(>\sqrt{n}\)的數\(x\)，\(n/x\)一定\(<\sqrt{n}\)”，就可以記下所有\(m\)的編號了。

int n,val[N*2],id1[N],id2[N];
//主函數（預處理）
int sqrt_n=sqrt(n),tot=0;
for(int i=1,j;i<=n;i=j+1) {
    j=n/(n/i);
    int w=n/i;val[++tot]=w;
    if(w<=sqrt_n) id1[w]=tot;
    else id2[n/w]=tot;
}
//查詢某個m的編號
inline int get_id(int m) {
    if(m<=sqrt_n) return id1[m];
    else return id2[n/m];
}

這一部分的時間復雜度被證明是\(O(\frac{n^{\frac{3}{4}}}{\log n})\)的。

STEP2 貢獻加和

現在我們用剛剛搞出來的\(G(m)=\sum_{i=1}^m[i\in prime]f(i)\)，求出所有\(f\)的和。

設\(S(n,j)=\sum_{i=1}^n[minp(i)\geq p[j]]f(i)\)。則顯然最終答案是\(f(1)+S(n,1)\)。

正如一開始所說的，我們分別考慮質數和合數對\(S(n,j)\)的貢獻。顯然質數的貢獻是\(G(n)-\sum_{i=1}^{j-1}f(p[i])\)。對於合數的貢獻，我們可以枚舉這個合數的“最小質因子及其次數”，然后進行遞推。容易列出遞推式：

\[S(n,j)=G(n)-\sum_{i=1}^{j-1}f(p[i])+\sum_{i=j}^{p[i]^2\leq n}\sum_{e=1}^{p[i]^{e+1}\leq n}\Big(f(p[i]^e)S(\lfloor\frac{n}{p[i]^e}\rfloor,i+1)+f(p[i]^{e+1})\Big) \]

式子里的二重循環就是在枚舉合數的“最小質因子及其次數”，注意合數可以有多個不同的質因子（繼續遞歸）或只有一個質因子的若干次方，這些都在式子里得以體現，請讀者仔細理解。

在求\(S\)時我們直接遞歸，不用記憶化。這一部分的時間復雜度，被證明也是\(O(\frac{n^{\frac{3}{4}}}{\log n})\)。反正你只要知道它能跑很大的數據（\(10^{10}\)的數據大約1秒搞定）就可以了。

一些栗子

wqy大神的良心模板題

預處理\(f_1(i)=i^2\)和\(f_2(i)=i\)在質數處的前綴和，在需要用到時減一下即可。

參考代碼：

//P5325
#include <bits/stdc++.h>
#define ll long long
#define pb push_back
#define mk make_pair
#define pii pair<int,int>
#define fst first
#define scd second
using namespace std;
inline ll read() {
	ll f=1,x=0;char ch=getchar();
	while(!isdigit(ch)){if(ch=='-')f=-1;ch=getchar();}
	while(isdigit(ch)){x=x*10+ch-'0';ch=getchar();}
	return x*f;
}
inline void write(ll x) {
	if(!x)putchar('0');if(x<0)x=-x,putchar('-');
	static int sta[20];register int tot=0;
	while(x)sta[tot++]=x%10,x/=10;
	while(tot)putchar(sta[--tot]+48);
}
const int MAXN=1e6+5;
const int MOD=1e9+7,INV6=166666668;
ll n,sqrtn,prm[MAXN],sump[MAXN],sum2p[MAXN],cnt,g1[MAXN],g2[MAXN];
bool v[MAXN];
void sieve(int n) {
	cnt=0;
	for(int i=2;i<=n;++i) {
		if(!v[i]) {
			prm[++cnt]=i;
			sump[cnt]=(sump[cnt-1]+i)%MOD;
			sum2p[cnt]=(sum2p[cnt-1]+(ll)i*i%MOD)%MOD;
		}
		for(int j=1;j<=cnt && prm[j]*i<=n;++j) {
			v[prm[j]*i]=1;
			if(i%prm[j]==0) break;
		}
	}
}
ll val[MAXN];
int id1[MAXN],id2[MAXN];
inline int get_id(ll x){return ((x<=sqrtn)?id1[x]:id2[n/x]);}
ll getS(ll x,ll y) {
	if(x<prm[y] || x<=1) return 0;
	int k=get_id(x);
	ll res=((g2[k]-g1[k]+MOD)%MOD-(sum2p[y-1]-sump[y-1]+MOD)%MOD+MOD)%MOD;
	for(int i=y;i<=cnt && prm[i]*prm[i]<=x;++i) {
		ll t1=prm[i],t2=prm[i]*prm[i];
		for(int j=1;t2<=x;++j,t1=t2,t2*=prm[i]) {
			ll tt1=t1%MOD,tt2=t2%MOD;
			res=(res+getS(x/t1,i+1)*tt1%MOD*(tt1-1)%MOD+tt2*(tt2-1)%MOD)%MOD;
		}
	}
	return res%MOD;
}
int main() {
	ios::sync_with_stdio(0);/*syn加速*/
	n=read();
	sieve(sqrtn=sqrt(n));
	int tot=0;/*O(sqrt(n))級別*/
	for(ll i=1,j;i<=n;i=j+1) {
		j=n/(n/i);
		ll w=n/i;
		val[++tot]=w;
		if(w<=sqrtn) id1[w]=tot;
		else id2[n/w]=tot;
		
		w%=MOD;
		g1[tot]=w*(w+1)/2%MOD;
		g2[tot]=w*(w+1)%MOD*(2LL*w+1)%MOD*INV6%MOD;
		
		g1[tot]=(g1[tot]-1+MOD)%MOD;
		g2[tot]=(g2[tot]-1+MOD)%MOD;
	}
	for(int j=1;j<=cnt;++j) {
		for(int i=1;i<=tot && prm[j]*prm[j]<=val[i];++i) {
			int k=get_id(val[i]/prm[j]);
			g1[i]=(g1[i]-prm[j]*(g1[k]-sump[j-1]+MOD)%MOD+MOD)%MOD;
			g2[i]=(g2[i]-prm[j]*prm[j]%MOD*(g2[k]-sum2p[j-1]+MOD)%MOD+MOD)%MOD;
		}
	}
	//for(int i=1;i<=tot;++i)cout<<val[i]<<" "<<g1[i]<<" "<<g2[i]<<endl;
	write((getS(n,1)+1)%MOD);puts("");
	return 0;
}

意外搞定了杜教篩的模板

先看\(\sum_{i=1}^n\varphi(i)\)。首先\(\varphi(p)=p-1\quad (p\in prime)\)。於是我們預處理\(f_1(i)=i\)和\(f_2(i)=1\)的前綴和，按與上一題相同的方法相減即可求出\(G(m)\)。

另外，\(\varphi(p^c)=p^{c-1}(p-1)\)（自己yy可得），這樣在求\(S\)時，質數的整次冪的函數值也很好求。

然后是\(\mu(i)\)，顯然\(\mu(p)=-1\quad (p\in prime)\)。於是我們只要對前面已經求好的\(f_2(i)=1\)求個相反數即可。

另外，\(\mu(p^k)=0\quad (k>1)\)，因此求\(S\)時不需要枚舉最小質因子的次數。

參考代碼：

//P4213 min_25 
#include <bits/stdc++.h>
#define ll long long
#define pb push_back
#define mk make_pair
#define pii pair<int,int>
#define fst first
#define scd second
using namespace std;
/*  -----  by:duyi  -----  */
const int N=5e4;
const int MAXN=N*2+5;
int n,sn,cnt,tot,p[MAXN],sum[MAXN],id1[MAXN],id2[MAXN],val[MAXN];
ll g1[MAXN],g2[MAXN];
bool v[MAXN];
void sieve() {
	v[1]=1;
	for(int i=2;i<=N;++i) {
		if(!v[i]) p[++cnt]=i;
		for(int j=1;j<=cnt && (ll)i*p[j]<=N;++j) {
			v[i*p[j]]=1;
			if(i%p[j]==0) break;
		}
	}
	for(int i=1;i<=cnt;++i) sum[i]=sum[i-1]+p[i];
}
inline int get_id(int x) {
	if(x<=sn) return id1[x];
	else return id2[n/x];
}
ll S_phi(int x,int y) {
	if(x<=1 || p[y]>x) return 0;
	ll res=g1[get_id(x)]-g2[get_id(x)]-(sum[y-1]-(y-1));
	for(int i=y;i<=cnt && (ll)p[i]*p[i]<=x;++i) {
		ll pre=1,cur=p[i];
		for(int j=1;cur*p[i]<=x;++j) {
			res+=pre*(p[i]-1LL)*S_phi(x/cur,i+1)+cur*(p[i]-1LL);
			pre=cur;cur=cur*p[i];
		}
	}
	return res;
}
ll S_mu(int x,int y) {
	if(x<=1 || p[y]>x) return 0;
	ll res=-g2[get_id(x)]+y-1;
	for(int i=y;i<=cnt && (ll)p[i]*p[i]<=x;++i) {
		res+=(-S_mu(x/p[i],i+1));
	}
	return res;
}
void solve() {
	cin>>n;
	sn=sqrt(n);tot=0;
	for(int i=1,j;i<=n;i=j+1) {
		j=n/(n/i);int w=n/i;
		val[++tot]=w;
		if(w<=sn) id1[w]=tot;
		else id2[n/w]=tot;
		g1[tot]=(ll)w*(w+1LL)/2LL-1LL;
		g2[tot]=w-1;
	}
	for(int i=1;i<=cnt;++i) {
		for(int j=1;j<=tot && (ll)p[i]*p[i]<=val[j];++j) {
			int t=get_id(val[j]/p[i]);
			g1[j]-=(ll)p[i]*(g1[t]-sum[i-1]);
			g2[j]-=(g2[t]-(i-1));
		}
	}
	cout<<S_phi(n,1)+1LL<<" "<<S_mu(n,1)+1LL<<endl;
}
int main() {
	ios::sync_with_stdio(0);/*syn加速*/
	sieve();
	int T;cin>>T;while(T--)solve();
	return 0;
}

總結一下

用min_25做題主要是要想清楚兩個問題：

1. 該函數在質數處的取值怎么求和？（找到合適的\(f'\)）

2. 質數的次冪處的取值\(f(p^c)\)能否快速求？

本文寫得比較粗淺，歡迎大家補充。如有紕漏歡迎指正。