min_25筛（学习笔记）

最新推荐文章于 2024-07-14 21:17:23 发布

QAQQQQQQQQQQQ

最新推荐文章于 2024-07-14 21:17:23 发布

阅读量2.1k

点赞数 3

分类专栏： min_25筛笔记文章标签： min_25筛

本文链接：https://blog.csdn.net/sizeof_you/article/details/86516018

版权

笔记同时被 2 个专栏收录

11 篇文章 0 订阅

订阅专栏

min_25筛

1 篇文章 0 订阅

订阅专栏

之前考试就学了一下 $min\_25$ 筛，现在理解得更多一些，补一波笔记

首先要知道这是一种可以筛出积性函数前缀和的东西，大概可以做洲阁筛能做的东西，复杂度是 $O(\frac{n^{\frac{3}{4}}}{logn})$ 的，大概能做到 $10^{10}$ ~ $10^{11}$

算法步骤

设 $f (i)$ 是一个积性函数，或者说 $f (i)$ 是一个关于质数 $p$ 的多项式，并且 $f(p^q)$ 可以快速求出来
要求 $\sum_{i=1}^nf(i)$

（下面内容中 $minp_x$ 表示 $x$ 的最小质因子， $p_j$ 表示筛出来的第 $j$ 个质数

线性筛出 $\sqrt{n}$ 以内的所有质数和 $f (i)$
求出 $\sum_{i=1}^n[i\in prime]f(i)$
设 $g (n, j)$ 表示 $\sum_{i=1}^n[i\in pirme\ or\ minp_i>p_j]f(i)$
考虑从 $g (n, j - 1) - > g (n, j)$ ，需要把 $p_j$ 的贡献算出来减掉，就是
$g(n,j)=\begin{cases}g(n,j-1)&p_j^2>n\\g(n,j-1)-f(p_j)\times [g(\frac{n}{p_j},j-1)-\sum_{i=1}^{j-1}f(p_i)]&p_j^2\le n\end{cases}$
上式就是考虑 $minp=p_j$ 的那些数的贡献，因为是积性函数，所以可以把 $p_j$ 提出来，剩下有
$g(\frac{n}{p_j},j-1)$ 那么多，其中不包括质数的所以要减掉 $\sum_{i=1}^{j-1}f(p_i)$
求出 $\sum_{i=1}^nf(i)$
这一步中其实只需要考虑合数的情况，因为质数已经都算出来了
设 $S(n,j)=\sum_{i=1}^n[minp_i\ge p_j]f(i)$ ，考虑每个合数在 $minp>p_j$ 时的贡献，可以枚举 $m i n p$ 和指数 $e$ ，就有
$S(n,j)=g(n,j)-\sum_{i=1}^{j-1}f(p_i)+\sum_{k=j}^{p_k^2\le n}\sum_{e=1}^{p_k^{e+1}\le n}S(\frac{n}{p_k^e},k+1)\times f(p_k^e)+f(p_k^{e+1})$
答案就是 $S(n,1)+f(1)=\sum_{i=1}^nf(i)$
这个通常是递归求解，也有非递归版的，有时会比递归版快很多，下面有

例题

LOJ#6053. 简单的函数
大概是 $min\_25$ 筛模板题？
发现除了 $p = 2$ 是 $f (p) = 3$ ，其余情况 $f (p) = p - 1$ ，所以算出 $g(n)=\sum_{i=1}^n[i\in prime]i$ 和 $h(n)=\sum_{i=1}^n[i\in prime] 1$ 就可以算出所有 $f (p)$ 的前缀和，进而算出 $S (n, 1)$ ，另外要特判一下含 $2$ 的情况，要 $+ 2$

代码如下：

#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<cmath>
#define N 1000005
#define LL long long
using namespace std;

const int mod=1e9+7;
LL n,w[N];
int cnt,pri[N],sqr,tot,m,id1[N],id2[N],g[N],h[N],sum[N];
bool vis[N];

inline int qpow(int x,int k){
	int ret=1;
	while(k){
		if(k&1) ret=1LL*ret*x%mod;
		x=1LL*x*x%mod; k>>=1;
	} return ret;
}

inline void get_prime(int n){
	for(int i=2;i<=n;i++){
		if(!vis[i]) pri[++cnt]=i,sum[cnt]=(sum[cnt-1]+i)%mod;
		for(int j=1;j<=cnt && i*pri[j]<=n;j++){
			vis[i*pri[j]]=1;
			if(i%pri[j]==0) break;
		}
	}
}

LL S(LL x,int y){
	if(x<=1 || pri[y]>x) return 0;
	int pos=(x<=sqr)?id1[x]:id2[n/x];
	LL ret=g[pos]-h[pos]-sum[y-1]+(y-1); ret=(ret+mod)%mod;
	if(y==1) ret+=2;
	for(int i=y;i<=cnt && 1LL*pri[i]*pri[i]<=x;i++){
		LL p1=pri[i],p2=1LL*pri[i]*pri[i];
		for(int e=1;p2<=x;e++,p1=p2,p2*=pri[i])
			(ret+=1LL*S(x/p1,i+1)*(pri[i]^e)%mod+(pri[i]^(e+1))%mod)%=mod;
	} return ret;
}

int main(){
	scanf("%lld",&n); int tmp=qpow(2,mod-2);
	sqr=sqrt(n); get_prime(sqr);
	for(LL i=1,j;i<=n;i=j+1){
		j=n/(n/i); w[++m]=n/i;
		if(w[m]<=sqr) id1[w[m]]=m;
		else id2[n/w[m]]=m;
		h[m]=(w[m]-1)%mod;
		g[m]=((w[m]+2)%mod)*((w[m]-1)%mod)%mod*tmp%mod;
	}
	for(int j=1;j<=cnt;j++)
		for(int i=1;i<=m&&1LL*pri[j]*pri[j]<=w[i];i++){
			int k=(w[i]/pri[j]<=sqr)?id1[w[i]/pri[j]]:id2[n/(w[i]/pri[j])];
			g[i]=(g[i]-1LL*pri[j]*(g[k]-sum[j-1])%mod+mod)%mod;
			h[i]=(h[i]-(h[k]-(j-1))+mod)%mod;
		}
	printf("%lld\n",S(n,1)+1);
	return 0;
}

UOJ#188. 【UR #13】Sanrd
$f (x)$ 表示 $x$ 的次大质因子，看起来不是积性函数，并且也不是关于 $p$ 的多项式

但它也是可以用 $min\_25$ 筛的！！！

首先可以考虑上面说的 $S$ 函数怎么求，是要枚举 $minp=p_j,p_j^e$ 考虑这些数的贡献，那实际上这个问题也可以这么解决，只考虑质因数 $\le 2$ 的时候，因为 $> 2$ 可以继续递归求解，如果次小质因子 $p_j$ ，那么贡献其实就是剩下的数中是 $p r i m e$ 且 $\ge p_j$ 的个数，然后就是前面那种套路的想法了
$S(n,j)=\sum_{k=j}^{p_k^2\le n}\sum_{e=1}^{p_k^{e+1}\le n}S(\frac{n}{p_k^e},k+1)+p_k\times \sum_{i=p_k}^{\frac{n}{p_k^e}}[i\in prime]$

代码如下：

#include<iostream>
#include<cstdio>
#include<cstring>
#include<cmath>
#include<algorithm>
#define N 1000005
#define LL long long
using namespace std;

LL l,r,w[N],g[N],n;
int cnt,pri[N],m,id1[N],id2[N],sqr;
bool vis[N];

inline void get_prime(int n){
	for(int i=2;i<=n;i++){
		if(!vis[i]) pri[++cnt]=i;
		for(int j=1;j<=cnt && i*pri[j]<=n;j++){
			vis[i*pri[j]]=1;
			if(i%pri[j]==0) break;
		}
	}
}

inline LL S(LL x,int y){
	if(x<=2||pri[y]>x) return 0;
	LL ret=0;
	for(int i=y;i<=cnt && 1LL*pri[i]*pri[i]<=x;i++)
		for(LL t=pri[i];t*pri[i]<=x;t*=pri[i]){
			int pos=(x/t<=sqr)?id1[x/t]:id2[n/(x/t)];
			ret+=S(x/t,i+1)+1LL*pri[i]*(g[pos]-(i-1));
	}
	return ret;
}

inline LL solve(LL nn){
	sqr=sqrt(nn); m=0; n=nn;
	for(LL i=1,j;i<=n;i=j+1){
		j=n/(n/i); w[++m]=n/i;
		if(w[m]<=sqr) id1[w[m]]=m;
		else id2[n/w[m]]=m;
		g[m]=w[m]-1;
	}
	for(int i=1;i<=cnt && 1LL*pri[i]*pri[i]<=n;i++)
		for(int j=1;j<=m && 1LL*pri[i]*pri[i]<=w[j];j++){
			int pos=(w[j]/pri[i]<=sqr)?id1[w[j]/pri[i]]:id2[n/(w[j]/pri[i])];
			g[j]-=g[pos]-(i-1);
		}
	return S(n,1);
}

int main(){
	scanf("%lld%lld",&l,&r);
	get_prime((int)sqrt(r));
	printf("%lld\n",solve(r)-solve(l-1));
	return 0;
}

LOJ#572. 「LibreOJ Round #11」Misaka Network 与求和

这道题其实就是上一道题套了个莫比乌斯反演
套路反演得到一个可以数论分块的式子，然后中间那部分再两边卷一下 $1$ 函数， $\mu*1=e$ ，右边就没了，然后左边移项就可以得到
$g(n)=\sum_{i=1}^nF(i)-\sum_{i=2}^ng(\frac{n}{i})$
发现 $F$ 前缀和就用上一个题的方法 $min\_25$ 筛筛出来，然后再数论分块一下就可以算这个式子了

递归版复杂度玄学，但可以跑过
非递归版复杂度十分优秀，大概就是 $O(min\_25筛)$ 级别的

代码如下：
（中间注释掉那部分是递归版的，非递归好像能跑递归的 $\frac{1}{3}$

#include<bits/stdc++.h>
#define N 1000005
#define uint unsigned int
using namespace std;

int K,m,cnt,pri[N],id1[N],id2[N],sqr;
bool vis[N];
uint n,w[N],g[N],f[N],h[N],ans,sum[N];

template<class T>inline void rd(T &x){
	x=0; short f=1; char c=getchar();
	while(c<'0' || c>'9') f=c=='-'?-1:1,c=getchar();
	while(c<='9' && c>='0') x=x*10+c-'0',c=getchar();
	x*=f;
}

inline int get(uint x){return x<=sqr?id1[x]:id2[n/x];}

inline uint qpow(uint x,int k){
	uint ret=1;
	while(k){
		if(k&1) ret*=x;
		x*=x; k>>=1;
	} return ret;
}

inline void get_prime(int n){
	for(int i=2;i<=n;i++){
		if(!vis[i]) pri[++cnt]=i,f[cnt]=qpow(i,K);
		for(int j=1;j<=cnt&&i*pri[j]<=n;j++){
			vis[i*pri[j]]=1;
			if(i%pri[j]==0) break;
		}
	}
}

inline void init(){
	for(uint i=1,j;i<=n;i=j+1){
		j=n/(n/i); w[++m]=n/i;
		if(w[m]<=sqr) id1[w[m]]=m;
		else id2[n/w[m]]=m;
		g[m]=w[m]-1;
	}
	for(int i=1;i<=cnt && pri[i]<=sqr;i++)
		for(int j=1;j<=m&&1LL*pri[i]*pri[i]<=w[j];j++)
			g[j]-=g[get(w[j]/pri[i])]-(i-1);
	
	for(int i=cnt;i;i--)
		for(int j=1;j<=m && 1LL*pri[i]*pri[i]<=w[j];j++)
			for(uint t=pri[i];1LL*t*pri[i]<=w[j];t*=pri[i])
				sum[j]+=sum[get(w[j]/t)]+f[i]*(g[get(w[j]/t)]-i+1);
}
/*
inline uint S(uint x,int y){
	if(x<=1 || pri[y]>x) return 0;
	uint ret=0;
	for(int i=y;i<=cnt && 1LL*pri[i]*pri[i]<=x;i++){
		uint p1=pri[i],p2=1LL*pri[i]*pri[i];
		for(uint t=pri[i];t*pri[i]<=x;t*=pri[i])
			ret+=S(x/t,i+1)+f[i]*(g[get(x/t)]-(i-1));
	} return ret;
}
*/
inline uint solve(uint x){
	int pos=get(x);
	if(h[pos]!=-1) return h[pos];
	uint ret=sum[pos]+g[pos];//minp>=pri[1] & prime
	for(uint i=2,j;i<=x;i=j+1){
		j=x/(x/i); ret-=solve(x/i)*(j-i+1);
	} return h[pos]=ret;
}

int main(){
	rd(n); rd(K); sqr=sqrt(n);
	memset(h,-1,sizeof h);
	get_prime(sqr); init();
	uint lst=0,cur,x;
	for(uint i=1,j;i<=n;i=j+1){
		x=n/i,j=n/x; cur=solve(j);
		ans+=(cur-lst)*x*x; lst=cur;
	}
	printf("%u\n",ans);
	return 0;
}