首页 > 技术文章 > 浅谈 Min_25 筛

Appleblue17 2020-12-24 18:15 原文

2022/7/5 Upd:略微修改了版面。

简介

和杜教筛一样,\(\text{Min\_25}\) 筛能够快速计算(快过线性)某些积性函数的前缀和。

\(\text{Min\_25}\) 筛的特点是最小质因数

跟这个有关的比较适合拿 \(\text{Min\_25}\) 筛解决。

最简单的形式是:\(f(p^e)\) 可以表示成关于 \(p^e\) 的低阶多项式。

接下来以这道题为例:【模板】Min_25筛


算法

题意

积性函数 \(f(x)\) 满足 \(f(p^k)=p^k(p^k-1)\) ,求:

\[\sum\limits_{i=1}^n f(i) \]

\(n \leq 10^{10}\)


首先定义几个东西

  • \(p_1,p_2,\dots,p_n\) 表示 前 \(n\) 个素数。

  • \(\mathbb P\) 表示 素数集。

  • \(low(n)\) 表示 \(n\) 的最小质因数为 \(p_{low(n)}\)

  • \(t\) 表示不超过 \(\sqrt{n}\) 的质数数量。


S(m,k)

  • 先提一点,由于 \(1\) 太特殊了,在计算过程中都不考虑 \(1\) ,最后再在答案里加上。

首先就是 \(\text{Min\_25}\) 筛最核心的部分:它是采用类似 DP 的方法求解的。

\(S(m,k)\) 为:不超过 \(m\) ,且最小质因数大于 \(p_k\) 的数的函数值之和。

或者说: \(S(m,k)=\sum\limits_{i=2}^m f(i)[low(i)>k]\)

那我们的答案就是 \(ans=S(n,0)\) 了。

好的,那怎么转移呢?

非常经典的转移方式:枚举最小质因数及它的次数。

比如说最小质因数为 \(p_j\) ,次数为 \(e\)

那我们就知道把 \(p_j^e\) 提出来,别的质因数一定要比 \(j\) 大了,也就是 \(S(\lfloor\frac{m}{p_j^e}\rfloor,j)\)

于是就可以转移啦……?

好,来试试!

\[S(m,k)=\sum\limits_{k<j} \sum\limits_{e \geq 1} f(p_j^e)(S(\lfloor\frac{m}{p_j^e}\rfloor,j)+1) \]

  • 解释:\(f(p_j^e)\) 是它的贡献(由于积性函数可以拎出来);\(+1\) 是指它自身(什么都不乘)有一个贡献。

成功啦。


等一下,这东西哪来的复杂度保证啊?

枚举所有质数能快就怪了。

咋办?

仔细一看,很容易发现有个限制:\(p_j^e \leq m\)

也就是说, \(p_j\) 比较大的时候,\(e=1\)

那我们把 \(e\) 只能为 \(1\) 的拎出来不就快多了吗?

先要算 \(p_k+1\)\(m\)质数的贡献。

于是式子可以改写一下:

\[S(m,k)=\sum\limits_{i=p_k+1}^m[i \in \mathbb P]f(i)+\sum\limits_{k<j \leq t,p_j^2 \leq m} \sum\limits_{e \geq 1}f(p_j^e)(S(\lfloor \frac{m}{p_j^e}\rfloor,j)+[e \neq 1]) \]

解释:当 \(p_j^2>m\) 时,\(S(\lfloor\frac{m}{p_j^e}\rfloor,j)=0\) 没有贡献,所以不用枚举。

好,接下来处理前面那个东西吧!

g(m,k)

这显然又可以改一下:

\[\sum\limits_{i=p_k+1}^m[i \in \mathbb P]f(i)=\sum\limits_{i=1}^m[i \in \mathbb P]f(i)-\sum\limits_{i=1}^{p_k}[i \in \mathbb P]f(i) \]

后面那个东西就是前 \(k\) 个质数的贡献,线性筛直接就得到了。

为了方便把它记成 \(s_k=\sum\limits_{i=1}^{k} f(p_i)\)

前面那个……挺不好求的。

为了方便把它记成 \(g_k=\sum\limits_{i=1}^{k}[i \in \mathbb P] f(i)\)


这时候就要请出第二个 DP:

\(g(m,k)\) 为:不超过 \(m\) ,且 最小质因数大于 \(p_k\) 的数 或者 质数 的函数值之和。

或者说: \(g(m,k)=\sum\limits_{i=2}^m f(i)[(low(i)>k) \vee (i \in \mathbb P)]\)

这跟 \(S(m,k)\) 挺相似的:\(g(m,k)=S(m,k)+s_k\)

那要求的就是 \(g_k=g(k,t)\)

转移?

这次我们换个转移方式,考虑从 \(g(m,k-1)\) 转移到 \(g(m,k)\)

显然少掉的只有 最小质因数为 \(p_k\) 的合数

(注意:上面 \(S\) 的计算中不这样转移是因为难以去掉大质数来优化)

啊,更简单了,连最小质因数都不用枚举。

不过呢,也不需要枚举次数 \(e\)

之前不是提到 \(f(x)\) 可以表示成低阶多项式吗?

拆开不就成 完全积性函数 了?

可以这么拆是因为现在只需要求 \(p^k\) 处的值

好,那假装现在 \(f(x)\) 变成了 \(x^k\)

\[g(m,k)=g(m,k-1)-f(p_k)(g(\lfloor\frac{m}{p_k}\rfloor,k-1)-g(p_{k-1},k-1)) \]

当然后面那个东西就是 \(s_{k-1}\)

\[g(m,k)=g(m,k-1)-f(p_k)(g(\lfloor\frac{m}{p_k}\rfloor,k-1)-s_{k-1}) \]

计算

注意到 \(g(m,k)\) 只会由 \(g(?,k-1)\) 转移而来。

而我们最后要用的只有 \(g(?,t)\),也许可以预处理。

再看回 \(S(m,k)\)

\[S(m,k)=g_m-s_k+\sum\limits_{k<j\leq t} \sum\limits_{e}f(p_j^e)(S(\lfloor\frac{m}{p_j^e}\rfloor,j)+[e \neq 1]) \]

我们只求 \(S(n,0)\)

然后你仔细看这个 \(m\):它从 \(n\) 开始,被整除了一个数,又被整除了一个数……

所以 \(m\) 必可以表示成 \(\lfloor\frac{n}{d}\rfloor\)\(d\) 是任意数)

这个只有 \(2\sqrt{n}\) 个,不错。

这是非常关键的一步,原本无法处理的问题通过这个根号降到了可以处理的范围。

(顺带一提,这个思想在杜教筛其实也是可以用的)


所以 \(g(?,t)\) 也只需要预处理出来 \(g(\lfloor\frac{n}{d}\rfloor,t)\) 即可。

当然这个预处理就是个 DP 啦。

然后把 \(\lfloor\frac{n}{d}\rfloor\) 离散化一下,套套式子:

\[g(m,k)=g(m,k-1)-f(p_k)(g(\lfloor\frac{m}{p_k}\rfloor,k-1)-s_{k-1}) \]

显然 \(k\) 这一维可以被滚动掉,就很优了。

当然注意 \(p_k^2 \leq m\)

求解直接递归即可,不用记忆化。

\[S(m,k)=g_m-s_k+\sum\limits_{k<j\leq t} \sum\limits_{e}f(p_j^e)(S(\lfloor\frac{m}{p_j^e}\rfloor,j)+[e \neq 1]) \]


再来介绍一个离散化的小技巧。

对于 \(k=\lfloor\frac{n}{d}\rfloor\) ,如果它不超过 \(\sqrt{n}\) ,就扔进 \(ind1[k]\) 里;否则就扔进 \(ind2[\lfloor\frac{n}{k}\rfloor]\)

好了,恭喜你学会了 \(\text{Min\_25}\) 筛!

代码

#include<bits/stdc++.h>
using namespace std;
const long long N=2e5+5,mod=1e9+7;
long long n,m;
long long pri[N],cnt;
long long mp[N],g1[N],g2[N],id;
long long s1[N],s2[N],ind1[N],ind2[N];
bool pd[N];

long long getpos(long long x){
	if(x<=m) return ind1[x];
	else return ind2[n/x];
}

void init(long long lim){
	pd[1]=1;
	for(long long i=2;i<=lim;i++){
		if(!pd[i]){
			pri[++cnt]=i;
			s1[cnt]=(s1[cnt-1]+i)%mod;
			s2[cnt]=(s2[cnt-1]+i*i%mod)%mod;
		}
		for(long long j=1;j<=cnt && i*pri[j]<=lim;j++){
			pd[i*pri[j]]=1;
			if(i%pri[j]==0) break;
		}
	}
}

long long S(long long x,long long k){
	if(pri[k]>=mp[x]) return 0;
	long long tot=(g2[x]-s2[k]-(g1[x]-s1[k])+mod+mod)%mod;
	for(long long j=k+1;j<=cnt && pri[j]*pri[j]<=mp[x];j++){
		long long tmp=pri[j];
		for(long long e=1;tmp<=mp[x];e++,tmp*=pri[j]){
			long long tmpp=tmp%mod;
			tot=(tot+tmpp*(tmpp-1)%mod*(S(getpos(mp[x]/tmp),j)+(e!=1))%mod)%mod;
		}
	}
	return tot;
}

int main(){
	cin>>n;
	m=sqrt(n);
	init(m);
	
	for(long long l=1,r;l<=n;l=r+1){
		r=n/(n/l);
		long long tmp=n/l;
		mp[++id]=tmp;//离散化对应的数 
		if(tmp<=m) ind1[tmp]=id;//离散化对应的下标 
		else ind2[n/tmp]=id;
		tmp%=mod;
		g1[id]=(tmp*(tmp+1)/2%mod+mod-1)%mod;
		g2[id]=(tmp*(tmp+1)/2%mod*((tmp*2+1)%mod)%mod*333333336%mod+mod-1)%mod;//此处是g(i,0)的值 
	}
	
	for(long long k=1;k<=cnt;k++){
		for(long long i=1;i<=id && pri[k]*pri[k]<=mp[i];i++){//g(i,k)
			long long tmp=getpos(mp[i]/pri[k]);
			g1[i]=(g1[i]+mod-pri[k]*((g1[tmp]+mod-s1[k-1])%mod)%mod)%mod;
			g2[i]=(g2[i]+mod-pri[k]*pri[k]%mod*((g2[tmp]+mod-s2[k-1])%mod)%mod)%mod;
		}
	}//dp预处理
	
	cout<<(S(getpos(n),0)+1)%mod;
}

拓展

还是这个东西

\[S(m,k)=g_m-s_k+\sum\limits_{k<j\leq t} \sum\limits_{e}f(p_j^e)(S(\lfloor\frac{m}{p_j^e}\rfloor,j)+[e \neq 1]) \]

刚才用递归,只求出了 \(S(n,0)\)

如果想求出所有 \(S(\lfloor\frac{n}{d}\rfloor,0)\) 呢?

众所周知搜索可以写成 DP。


转移是 \(k\) 从大往小转移的。

考虑从 \(S(m,k)\) 转移到 \(S(m,k-1)\)

\[S(m,k-1)=g_m-s_{k-1}+\sum\limits_{k-1<j\leq t} \sum\limits_{e}f(p_j^e)(S(\lfloor\frac{m}{p_j^e}\rfloor,j)+[e \neq 1]) \]

\[S(m,k)=g_m-s_k+\sum\limits_{k<j\leq t} \sum\limits_{e}f(p_j^e)(S(\lfloor\frac{m}{p_j^e}\rfloor,j)+[e \neq 1]) \]

\[S(m,k-1)-S(m,k)=s_k-s_{k-1}+\sum\limits_{e}f(p_k^e)(S(\lfloor\frac{m}{p_k^e}\rfloor,k)+[e \neq 1]) \]

\[S(m,k-1)=S(m,k)+s_k-s_{k-1}+\sum\limits_{e}f(p_k^e)(S(\lfloor\frac{m}{p_k^e}\rfloor,k)+[e \neq 1]) \]

挺像 背包 的不是吗?

for(long long i=1;mp[i]>=pri[cnt];i++){
	ans[i]=((g2[i]-g1[i])-(s2[cnt]-s1[cnt])+mod+mod)%mod;
}
for(long long j=cnt;j>=1;j--){
	for(long long i=1;mp[i]>=pri[j];i++){
		ans[i]=(ans[i]+(s2[j]-s1[j])-(s2[j-1]-s1[j-1])+mod+mod)%mod;
		if(pri[j]*pri[j]>mp[i]) continue;
		long long tmp=pri[j];
		for(long long e=1;tmp<=mp[i];e++,tmp*=pri[j]){
			long long tmpp=tmp%mod;
			ans[i]=(ans[i]+tmpp*(tmpp-1)%mod*(ans[getpos(mp[i]/tmp)]+(e!=1)))%mod;
		}
	}
}

不过这么写还是挺慢了(所以没啥实际意义)。

但是如果 \(f\) 只有 \(1\) 项的话,可以用完全背包优化。

练习

LOJ6053 简单的函数

题意:积性函数 \(f(n)\) 满足 \(f(p^k)=p \oplus k\)\(f(1)=1\)

求:

\[\sum\limits_{i=1}^n f(i) \]

\(n \leq 10^{10}\)


一看就是 Min_25 模板。

这个异或看起来很难搞,但是 \(S(m,k)\) 函数的计算是不需要表示成多项式的。

所以只用考虑 \(g(m,k)\) 计算时质数的取值。

显然:

  • \(f(2)=1\)

  • \(f(p)=p+1,p \neq 2\)

就变成模板啦~


UOJ188 【UR #13】Sanrd

题意:令 \(f(n)\)\(n\) 的次大质因数(重复的也算,比如 \(f(4)=2\) )。若 \(n\) 为质数,则规定 \(f(n)=0\)

求:

\[\sum\limits_{i=l}^r f(i) \]

\(l+r \leq 10^{11}\)


还是根据 Min_25 的经典套路,设出 \(S(m,k)\)

然后发现因为质数的时候没贡献,不用考虑了!

也就不用设 \(g\) 了,简单了许多。

分类讨论 \(p_j\) 是不是次大质因数,发现还要快速求质数个数。

求质数个数的话令模板那里的 \(f(x)=1\) 即可。


LOJ572 「LibreOJ Round #11」Misaka Network 与求和

题意:令 \(f(n)\)\(n\) 的次大质因数(重复的也算,比如 \(f(4)=2\) )。若 \(n\) 为质数,则规定 \(f(n)=1\)

求:

\[\sum\limits_{i=1}^n \sum\limits_{j=1}^n f(gcd(i,j))^k \]

\(n \leq 2 \times 10^{9}\)


莫反+杜教+Min_25 套娃题。

先令 \(f(x)=f(x)^k\),因为 Min_25 一样可以处理 \(f(x)^k\)

首先有个 gcd ,自然莫反把它搞出来。

然后发现需要快速求 \(g=f*\mu\) 的前缀和,而且只需要求 \(\lfloor\frac{n}{d}\rfloor\) 处的。

莫反一下,得到 \(f=g*1\)

杜教筛上,现在只要需要求 \(f\) 的前缀和了,而且还是只需要求 \(\lfloor\frac{n}{d}\rfloor\) 处的。

就是刚刚的 \(\text{Min\_25}\) 筛。

不过这里需要求所有位置的值,要用刚刚那个拓展。

解决~

推荐阅读