Min_25 筛 学习笔记

前言

很久之前就应该学的,但是因为太鸽了咕到了现在,而且估计以后也会忘,所以写一篇符合自己理解方式的学习笔记,估计对萌新的学习没有什么帮助,所以这边先放几个学习链接。

wucstdio

yyb

AThousandMoon

关于本文中一些符号的定义

  • (lpf(i))表示(i)的最小质因子。

  • 在不进行特殊说明的情况下,(p)视为质数。

Min_25 筛是干什么的?

Min_25 筛是用来求一个积性函数前缀和的东西。也就是求解这样的式子:(sum_{i=1}^{n} f(i))(其中(f)是积性函数)

可以使用 Min_25 筛的前提是什么?

要求积性函数(f)(f(p^k))的位置上能够快速地求出来,换一句话说,就是(f(p^k))是一个次数比较低的多项式。

怎么进行 Min_25 筛?

假设我们需要求(sum_{i=1}^n f(i))

那么,我们将(i)按照其是否为质数分类,原式则可以变为:(sum_{1leq p leq n} f(p) +sum_{i=1}^{n} [i\, is\, not\, a\, prime]f(i))

前面的先放一边,因为可以直接算出点值,考虑后面的合数部分怎么化。

因为(f)是积性函数,所以枚举质数,接下来可以把式子化成这样(定义(lpf(i))(i)的最小质因子):
(sum_{1leq p leq n} f(p) + sum_{p^e leq n} f(p^e) (sum_{i=1 & lpf(i)>p}^{left lfloor frac{n}{p^e} ight floor} f(i)))

但是,我们注意到如果是合数的话,那么它的最小质因子一定小于(sqrt{n})

所以式子就变成了这样:
(sum_{1leq p leq n} f(p) + sum_{p^e leq n& 1leq p leq sqrt{n} } f(p^e) (sum_{i=1 & lpf(i)>p}^{left lfloor frac{n}{p^e} ight floor} f(i)))

接下来我们要分两个步骤处理。

步骤 1

为了解决这个问题,我们设计一个函数(g(n,j))表示在(n)以内质数或者最小质因子大于第(j)个质数的数字的(h)函数值的和,当然,这个(h)函数并不是原来的(f)函数,只是一个在(p^k)处的取值和(f)相同或者相似的函数,因为我们的前提中(f(p^k))是一个关于(p^k)的低阶多项式,所以(h)一般会设为(h(i)=i^k)这种形式(注:可以设多个(h),最后加起来)。

下面用(p_j)表示第(j)个质数。

那么,用数学式子来表示的话,就是(g(n,j)=sum_{i=1}^{n} [i\, is\, a\, prime\, or\, lpf(i) > p_j] i^k)

我们考虑用 DP 的方法来解决这个问题。

考虑从(g(n,j-1))转移到(g(n,j))

很明显,需要减去一些数,因为有一些数已经不满足条件了,考虑这些不满足条件的数是哪些。

很明显,所有不满足条件的数就是最小质因子是(p_j)的合数。

那么,这些数我们提取除一个(p_j)出来,我们就可以写出这么一个式子来表示这些数对原来答案的贡献:(p_j^kcdot (g(frac{n}{p_j},j-1)- g(p_{j-1},j-1))),后面减去(g(p_{j-1},j-1))是因为我们要把原来的质数给去掉。

接下来我们就可以列出转移的式子了。

(g(n,j)=g(n,j-1)-p_j^kcdot (g(frac{n}{p_j},j-1)- g(p_{j-1},j-1)))

这个式子看起来计算是(O(n*|P|))(|P|)表示(n)以内的质数个数),但是我们发现(left lfloor frac{n}{x} ight floor)最多只有(O(sqrt{n}))钟取值,所以可以离散化,那么就可以快速地计算出(g(n,j))了,空间上的优化直接用滚动数组就可以了。

步骤 2

但是我们并不能通过(g)函数来找出我们所需要的答案。

所以我们仿照上面的思路,构造出一个(S(n,k))函数表示从(1)(n)最小质因子大于(p_k)(f)函数值的和,那么,如果我们能够快速地求出(S(n,0)),这一道题就解决了。

同样是分质数和合数来考虑,合数也依然是枚举最小质因子。

(S(n,k)=g(n,|P|)-g(p_{|P|},|P|)+sum_{p_j^eleq n & j>k}) (f(p_{j}^{e}) (S(left lfloor frac{n}{p_{j}^{e}} ight floor ,k)+[e eq 1]))

好的,这样(S)就可以递归算了,时间复杂度分析我不会,据说是(O(frac{n^{frac{3}{4}}}{log n}))

下面附上洛谷模板的代码:

#include <cmath>
#include <cstdio>
template<typename Elem>
void read(Elem &a){
	a=0;
	char c=getchar();
	while(c<'0'||c>'9'){
		c=getchar();
	}
	while(c>='0'&&c<='9'){
		a=(a<<1)+(a<<3)+(c^48);
		c=getchar();
	}
}
template<typename Elem>
void write(Elem a){
	if(a<10){
		putchar(a+'0');
		return;
	}
	write(a/10);
	putchar(a%10+'0');
}
const int Maxn=200000;
const int Mod=1000000007;
typedef long long ll;
const int Inv_2=(Mod+1)>>1;
const int Inv_3=(Mod+1)/3;
ll n;
int sqr;
bool np[Maxn+5];
int p[Maxn+5],p_len;
int sp_1[Maxn+5],sp_2[Maxn+5];
int g_1[Maxn+5],g_2[Maxn+5];
int pos_1[Maxn+5],pos_2[Maxn+5];
ll w[Maxn+5];
int tot;
void init(){
	np[0]=np[1]=1;
	for(int i=2;i<=Maxn;i++){
		if(!np[i]){
			p[++p_len]=i;
			sp_1[p_len]=(sp_1[p_len-1]+i)%Mod;
			sp_2[p_len]=(sp_2[p_len-1]+1ll*i*i)%Mod;
		}
		for(int j=1,x;j<=p_len&&(x=i*p[j])<=Maxn;j++){
			np[x]=1;
			if(i%p[j]==0){
				break;
			}
		}
	}
}
int find_sum(ll x,int y){
	if(p[y]>=x){
		return 0;
	}
	int k=x<=sqr?pos_1[x]:pos_2[n/x];
	int ans=((g_2[k]-g_1[k]+Mod)%Mod-(sp_2[y]-sp_1[y]+Mod)%Mod+Mod)%Mod;
	for(int i=y+1;i<=p_len&&1ll*p[i]*p[i]<=x;i++){
		ll p_e=p[i];
		for(int j=1;p_e<=x;j++,p_e*=p[i]){
			int p_e_m=p_e%Mod;
			ans=(ans+1ll*p_e_m*(p_e_m-1)%Mod*(find_sum(x/p_e,i)+(j!=1)))%Mod;
		}
	}
	return ans;
}
int main(){
	init();
	read(n);
	sqr=(int)sqrt(n+0.0);
	for(ll i=1,j;i<=n;i=j+1){
		j=n/(n/i);
		w[++tot]=n/i;
		g_1[tot]=w[tot]%Mod;
		g_2[tot]=1ll*g_1[tot]*(g_1[tot]+1)/2%Mod*(2*g_1[tot]+1)%Mod*Inv_3%Mod-1;
		if(g_2[tot]<0){
			g_2[tot]+=Mod;
		}
		g_1[tot]=1ll*g_1[tot]*(g_1[tot]+1)/2%Mod-1;
		if(g_1[tot]<0){
			g_1[tot]+=Mod;
		}
		if(n/i<=sqr){
			pos_1[n/i]=tot;
		}
		else{
			pos_2[n/(n/i)]=tot;
		}
	}
	for(int i=1;i<=p_len;i++){
		for(int j=1;j<=tot&&1ll*p[i]*p[i]<=w[j];j++){
			int k=w[j]/p[i]<=sqr?pos_1[w[j]/p[i]]:pos_2[n/(w[j]/p[i])];
			g_1[j]=(g_1[j]-1ll*p[i]*(g_1[k]-sp_1[i-1]+Mod)%Mod+Mod)%Mod;
			g_2[j]=(g_2[j]-1ll*p[i]*p[i]%Mod*(g_2[k]-sp_2[i-1]+Mod)%Mod+Mod)%Mod;
		}
	}
	write((find_sum(n,0)+1)%Mod);
	putchar('
');
	return 0;
}
原文地址:https://www.cnblogs.com/withhope/p/13258511.html