Min_25 筛学习笔记

Min_25 筛学习笔记

事实上我又学习了一个有点春的筛法。Min_25 筛用于求解积性函数的前缀和,即形如 g(n)=i=1nf(i) 形式的函数 g

众所周知,朴素筛法之所以无法做到低于线性是因为枚举了区间内的每一个数,那么我们想要做到低于线性,就必然需要做到通过一种方法将所有数字分成两类,用过一类的求解辅助另一类的求解,不难想到质数和合数。那么我们分开考虑如何求解质数和合数。对于 1 这个既不是质数也不是合数的数,我们选择将其的贡献最后单独统计,毕竟一个积性函数要么所有位置都是 0,要么 f(1)=1

质数求解

首先考虑如何求解质数的贡献。众所周知,幂函数是一种完全积性函数,假若当 x=pf(x)=xi,那么我们便可以将每一项分开求解。更具体的我们尝试统计 iPik 的值。

埃筛是一种复杂度略差于线性筛的筛法,但只表现为拥有较大的常数,在这里我们借鉴埃筛的思路。考虑埃筛利用每一个范围内的质数进行过滤,更具体的,对于每一个质数 p,其都将 p2modp=0 的数筛掉。那么我们考虑所有 n 的合数一定有一个 n 的质因子,因此我们只需要筛出 n 以下的素数即可。但是埃筛的时间复杂度依旧不可接受。考虑用 dp 优化这个过程,我们设 gi,j 表示 [1,i] 之间的数在经过 j 轮筛之后所剩余的数的贡献,那么我们需要知道的就是 gn,|Pn|。首先有初始值

gn,0=i=2nik

那么考虑因为质数 pj 只会筛掉 pj2 的数,那么我们可以推断出转移方程

gi,j=gi,j1(i<pj2)

接着考虑 ipj2 的时候,因为完全积性函数的性质,所以我们单独将被筛掉的合数的质因子 pj 提出来。那么剩下的部分应该是最小质因子 pj 的部分,因为如果最小质因子 <pj 一定会被之前的某个素数筛掉,我们可以用 gipj,j1 来表示,但是这里面包含了 [1,pj1] 的质数,不应该减去这部分,因此还需要将这部分贡献加回来,所以有转移方程

gi,j=gi,j1pjk(gipj,j1gpj1,j1)(ipj2)

不难发现可以用滚动数组优化掉第二维,但是时空复杂度依旧没有消下去。首先后半部分其实没有必要利用 g 数组来表示,反倒可以利用素数前缀和数组 s,而这个数组后面也会遇到。接着考虑可能被转移的位置 ipj,如果从搜索的角度看,就是 n,np1,np1p2,,有一个化简技巧就是 nab=nab。那么我们发现有用的下标其实只有 na,那么我们可以利用数论分块在 O(n) 的时间复杂度内处理出所有可能的值。

合数求解

在求解完素数之后我们就要酝酿着统计合数对答案的贡献了。因为给定的函数是积性函数的缘故,因此我们可以将合数表示为最小质因子的次方和另一个数相乘的结果,因为后半部分的最小质因子不能小于原来的最小质因子,因此我们令 S(n,b) 表示范围 [1,n] 内所有最小质因子 >pb 的数的贡献和,那么所求即为 S(n,0)。递归求解这个东西。当 npb 时,直接返回 0 即可。否则我们先统计所有素数的答案,然后进行最小质因子和其次幂的枚举,并将对应的答案贡献累加,更具体的,有转移方程

S(n,b)=gn,|Pn|i=1bf(pi)+i=b+1|Pn|1j,pijn(S(npij,i)+[j1]f(pij))

我们发现中间的求和可以用之前的质数前缀和快速求解,而剩余部分自行递归即可。整体的复杂度是 O(n34logn)

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]-p[j]*(g[id]-sum[j-1]))%P;//直接用前缀和代替 g 数组
		}
	}
	cout<<(SumF(n,0)+1+P)%P;//记得单独加上 1 的贡献
	return 0;
}
posted @   DycIsMyName  阅读(19)  评论(0编辑  收藏  举报
相关博文:
阅读排行:
· 震惊!C++程序真的从main开始吗?99%的程序员都答错了
· 【硬核科普】Trae如何「偷看」你的代码?零基础破解AI编程运行原理
· 单元测试从入门到精通
· 上周热点回顾(3.3-3.9)
· winform 绘制太阳,地球,月球 运作规律
点击右上角即可分享
微信分享提示