min_25 筛入门

什么是 min_25 筛

min_25 筛和洲阁筛、杜教筛一样，是一种低于线性的用于求积性函数前缀和的筛法。常用 min_25 筛的时间复杂度为 $O(\frac{n^{\frac34}}{\log n})$ ，而经过优化可以达到 $O(n^{\frac23})$ （但是常数巨大且一般用不着）。

前置知识

数论函数

数论函数：满足 $f:\mathbb N_+\rightarrow \mathbb C$ 的函数为数论函数。一般用到的是 $f:\mathbb N_+\rightarrow \mathbb Z$ 的函数。

积性函数：若数论函数 $f(x)$ 满足 $\forall a,b\in \mathbb N_+,(a,b)=1\Rightarrow f(ab)=f(a)f(b)$ ，则称 $f(x)$ 为积性函数。

完全积性函数：若数论函数 $f(x)$ 满足 $\forall a,b\in \mathbb N_+\Rightarrow f(ab)=f(a)f(b)$ ，则称 $f(x)$ 完全积性函数。

狄利克雷卷积：对于数论函数 $f(x)$ 和 $g(x)$ ，它们的狄利克雷卷积 $(f*g)(x)=\sum_{i|x}f(i)g(\frac x i)$ 。

常见的积性函数有： $\varphi(x),\mu(x),$ ~~题目定义的~~ $,\dots\dots$

常见的完全积性函数有： $id(x)=x; \epsilon(x)=[x=1]; I(x)=1; ......$

一些有趣的性质：

\begin{matrix} (1) & μ * I = ϵ \end{matrix}

$\mu * I = \epsilon\tag{1}\\$

\begin{matrix} (2) & μ * i d = φ \end{matrix}

$\mu * id = \varphi\tag{2}$

以下， $P$ 为质数集合， $p_i(i>0)$ 表示第 $i$ 个质数， $mp(n)(n>0)$ 表示 $n$ 的最小质因子。

埃拉托色尼筛

算法思想：对于质数 $p$ ，将 $p$ 的倍数全部筛掉。时间复杂度为 $O(n\log\log n)$

注意：对于质数 $p$ ，被筛掉数一定 $\ge p^2$ ；否则它一定会有 $< p$ 质因子。

欧拉筛

算法思想：对于每个数 $i$ ，筛掉那些 $mp\le mp(i)$ 的数。可以发现每个数只会在 $mp$ 处被筛掉，因此是线性筛。

min_25 筛

一般筛不够优秀的原因就是，它们枚举了数据范围内的每一个数，因此时间不会低于线性。

而低于线性的筛，就是利用算数基本原理，先计算质数的贡献，再推出合数的贡献。

min_25 筛的计算条件为：积性函数 $f(x)$ 在质数处可以被表示为简单多项式，并且对于质数， $f(p^c)$ 可以被高速算出。

计算质数贡献

我们考虑先计算质数的贡献。即：

\sum_{i = 1}^{n} [i \in P] f (i)

$\sum_{i=1}^{n} [i\in P] f(i)$

由于 $f$ 在质数处可以被表示为简单多项式，所以我们可以对于多项式的每一项分别计算。即只用考虑：

\sum_{i = 1}^{n} [i \in P] i^{k}

$\sum_{i=1}^n [i\in P] i^k$

考虑这样一个 DP ：

$g_{n,j}$ ：前 $n$ 个数进行 $j$ 轮埃筛之后的贡献和。

可以表达为：

g_{n, j} = \sum_{i = 2}^{n} [n \in P \lor m p (n) > p_{j}] i^{k}

$g_{n,j}=\sum_{i=2}^n [n\in P\or mp(n)> p_j]i^k$

设 $P'=[1,\lfloor\sqrt n\rfloor]\cap P$ ，即 $\sqrt n$ 范围内的质数。不难发现，由于每个合数 $m$ 必然有一个 $\le \sqrt m$ 的质因子，因此我们只需要将 $P'$ 中的质数全部筛一遍， $[1,n]$ 剩下的都是质数了。故 $g_{n,|P'|}$ 就是我们需要的质数的贡献。

初始状态为：

g_{n, 0} = \sum_{i = 2}^{n} i^{k}

$g_{n,0}=\sum_{i=2}^n i^k$

可以得到转移为：

g_{n, j} = {\begin{cases} g_{n, j - 1} & a < p_{b}^{2} \\ g_{n, j - 1} - p_{b}^{k} (g_{⌊ \frac{a}{p_{b}} ⌋, b - 1} - g_{p_{b - 1}, b - 1}) & a \geq p_{b}^{2} \end{cases}

$g_{n,j}=\begin{cases}g_{n,j-1}& a< p_b^2\\g_{n,j-1}-p_b^k(g_{\lfloor\frac a{p_b}\rfloor,b-1} - g_{p_{b-1},b-1})&a\ge p_b^2\end{cases}$

当 $a< p_b^2$ 的时候，不会筛掉任何数；否则考虑减掉的贡献。显然那些最小质因子为 $p_b$ 的数会被筛掉，即 $g_{\lfloor\frac a{p_b}\rfloor,b-1}$ ；但是这个东西里面有 $1\sim b-1$ 的质数，不应该减掉，因此还要再补上 $g_{p_{b-1},b-1}$ 。

这一部分可以用滚动数组优化。

计算总贡献

对于总的贡献，我们设这样一个函数 $S(a,b)$ ：

S (a, b) = \sum_{i = 2}^{a} [m p (i) \geq p_{b}] f (i)

$S(a,b)=\sum_{i=2}^a [mp(i)\ge p_b] f(i)$

请注意这里并没有质数的专门贡献，且里面的要求是 " 最小质因子不小于 $p_b$ " 。

其中质数的贡献已经算出来了。对于合数的情况，我们枚举它的最小质因子和其指数，可以得到：

S (n, j) = g (n, | P^{'} |) - \sum_{i = 1}^{j - 1} f (p_{i}) + \sum_{i = j}^{| P^{'} |} \sum_{1 \leq e, p_{i}^{e + 1} \leq a} (S (⌊ \frac{a}{p_{i}^{e}} ⌋, i + 1) + f (p_{i}^{e + 1}))

$S(n,j)=g(n,|P'|)-\sum_{i=1}^{j-1}f(p_i)+\sum_{i=j}^{|P'|}\sum_{1\le e, p_i^{e+1}\le a}\left(S(\lfloor\frac a{p_i^e}\rfloor,i+1)+f(p_i^{e+1})\right)$

其中 $\sum_{i=1}^{b-1}f(p_i)$ 可以预先筛出来。然后 $S(\lfloor\frac a{p_i^e}\rfloor,i+1)$ 可以递归下去继续算。

实现

首先需要对 $[1,\sqrt n]$ 里面的数进行一发线性筛求出质数。

根据整除分块理论，我们在代入 $S(n,1)$ 计算的时候，实际上 $a$ 的取值只有 $O(\sqrt n)$ 个。因此我们可以预处理这 $O(\sqrt n)$ 个取值的离散化后下标。设其中一种取值为 $x$ ，那么当 $x\le \sqrt n$ 的时候，我们直接将下标存下来；否则，由于 $\lfloor\frac n x\rfloor\le \sqrt n$ ，我们将下标存在 $\lfloor\frac n x\rfloor$ 里面。

例题

[LOJ]区间素数个数

这里不需要求 $S$ 的步骤，直接用第一步就可以了。

代码如下：

#include <cmath>
#include <cstdio>

typedef long long LL;

const int MAXS = 1e6 + 5;

template<typename _T>
void read( _T &x )
{
	x = 0;char s = getchar();int f = 1;
	while( s > '9' || s < '0' ){if( s == '-' ) f = -1; s = getchar();}
	while( s >= '0' && s <= '9' ){x = ( x << 3 ) + ( x << 1 ) + ( s - '0' ), s = getchar();}
	x *= f;
}

template<typename _T>
void write( _T x )
{
	if( x < 0 ){ putchar( '-' ); x = ( ~ x ) + 1; }
	if( 9 < x ){ write( x / 10 ); }
	putchar( x % 10 + '0' );
}

LL val[MAXS << 1], g[MAXS];
int id1[MAXS], id2[MAXS];
int prime[MAXS], pn;
LL N;
int s, tot;
bool isPrime[MAXS];

void EulerSieve( const int siz )
{
	isPrime[1] = true;
	for( int i = 2 ; i <= siz ; i ++ )
	{
		if( ! isPrime[i] ) prime[++ pn] = i;
		for( int j = 1 ; j <= pn && 1ll * i * prime[j] <= siz ; j ++ )
		{
			isPrime[i * prime[j]] = true;
			if( ! ( i % prime[j] ) ) break;
		}
	}
}

int getID( const LL x ) { return x <= s ? id1[x] : id2[N / x]; }

signed main()
{
	read( N );
	s = sqrt( N ), tot = 0;
	EulerSieve( s );
	for( LL l = 1, r, v ; l <= N ; l = r + 1 )
	{
		v = N / l, r = N / ( N / l );
		if( v <= s ) id1[v] = ++ tot; 
		else id2[N / v] = ++ tot;
		val[tot] = v, g[tot] = v - 1;
	}
	for( int j = 1 ; j <= pn ; j ++ )
		for( int i = 1 ; i <= tot && 1ll * prime[j] * prime[j] <= val[i] ; i ++ )
			g[i] -= g[getID( val[i] / prime[j] )] - ( j - 1 );
	write( g[getID( N )] ), putchar( '\n' );
	return 0;
}

[LG P4213]杜教筛

正所谓 " 树套树的题怎么能用树套树做呢？ " ，杜教筛的题怎么能用杜教筛？

考虑 min_25 ——事实上， min_25 最重要的地方就是求出 $g$ ，而后推 $S$ 其实就是板子的事情了。

当 $p$ 为质数时， $\varphi(p)=p-1, \mu(p)=-1$ ，因此我们只需要用 $g$ 求出素数个数和素数和即可。之后就是递归的事情了。

代码如下：

#include <cmath>
#include <cstdio>

typedef long long LL;

#define int LL

const int MAXS = 1e5 + 5;

template<typename _T>
void read( _T &x )
{
	x = 0;char s = getchar();int f = 1;
	while( s > '9' || s < '0' ){if( s == '-' ) f = -1; s = getchar();}
	while( s >= '0' && s <= '9' ){x = ( x << 3 ) + ( x << 1 ) + ( s - '0' ), s = getchar();}
	x *= f;
}

template<typename _T>
void write( _T x )
{
	if( x < 0 ){ putchar( '-' ); x = ( ~ x ) + 1; }
	if( 9 < x ){ write( x / 10 ); }
	putchar( x % 10 + '0' );
}

template<typename _T>
_T MAX( const _T a, const _T b )
{
	return a > b ? a : b;
}

LL gPhi[MAXS], gMu[MAXS], g1[MAXS], g2[MAXS];
LL ps[MAXS], ms[MAXS];
int val[MAXS], id1[MAXS], id2[MAXS];
int prime[MAXS], pn;
int N, s, cnt, tot;
bool isPrime[MAXS];

LL sqr( const LL x ) { return x * x; }
int getID( const int x ) { return x <= s ? id1[x] : id2[N / x]; }

void EulerSieve( const int siz )
{
	isPrime[1] = true;
	for( int i = 2 ; i <= siz ; i ++ )
	{
		if( ! isPrime[i] ) prime[++ pn] = i;
		for( int j = 1 ; j <= pn && 1ll * i * prime[j] <= siz ; j ++ )
		{
			isPrime[i * prime[j]] = true;
			if( ! ( i % prime[j] ) ) break;
		}
	}
	for( int i = 1 ; i <= pn ; i ++ ) ps[i] = ps[i - 1] + prime[i] - 1, ms[i] = ms[i - 1] - 1;
}

LL SPhi( const int a, const int b )
{
	if( a < prime[b] ) return 0;
	LL ret = gPhi[getID( a )] - ps[b - 1];
	if( b > tot ) return a <= prime[tot] ? 0 : ret;
	LL phi, p, tmp;
	for( int i = b ; i <= tot && 1ll * prime[i] * prime[i] <= a ; i ++ )
	{
		phi = prime[i] - 1, p = prime[i];
		for( int j = 1 ; p * prime[i] <= a ; j ++, p *= prime[i], phi *= prime[i] )
			ret += ( SPhi( a / p, i + 1 ) * phi + p * ( prime[i] - 1 ) );
	}
	return ret;
}

LL SMu( const int a, const int b )
{
	if( a < prime[b] ) return 0;
	LL ret = gMu[getID( a )] - ms[b - 1];
	if( b > tot ) return a <= prime[tot] ? 0 : ret;
	for( int i = b ; i <= tot && 1ll * prime[i] * prime[i] <= a ; i ++ ) ret -= SMu( a / prime[i], i + 1 );
	return ret;
}

signed main()
{
	int T;
	read( T );
	EulerSieve( 1e5 );
	while( T -- )
	{
		read( N );
		s = sqrt( N );
		for( tot = 1 ; prime[tot] <= s ; tot ++ );
		tot --, cnt = 0;
		for( int l = 1, r, v ; l <= N ; l = r + 1 )
		{
			r = N / ( v = N / l );
			if( v <= s ) id1[v] = ++ cnt;
			else id2[N / v] = ++ cnt;
			val[cnt] = v;
			g1[cnt] = v - 1, g2[cnt] = ( 1ll * v * ( v + 1 ) >> 1 ) - 1;
		}
		for( int j = 1, k ; j <= tot ; j ++ )
			for( int i = 1 ; i <= cnt && 1ll * prime[j] * prime[j] <= val[i] ; i ++ )
			{
				k = val[i] / prime[j];
				g1[i] -= g1[getID( k )] - ( j - 1 );
				g2[i] -= 1ll * prime[j] * ( g2[getID( k )] - g2[getID( prime[j - 1] )] );
			}
		for( int i = 1 ; i <= cnt ; i ++ ) gPhi[i] = g2[i] - g1[i], gMu[i] = - g1[i];
		write( SPhi( N, 1 ) + 1 ), putchar( ' ' ); 
		write( SMu( N, 1 ) + 1 ), putchar( '\n' );
	}
	return 0;
}

历史遗留问题

这是在 2022 的 update，距离我写这篇博客大约有两年。这段时间内基本上没有写过这个算法，~~主要还是因为 Powerful Numbers 这套工具太香了~~，所以冒出来了一些问题。

求 $g$ 和求 $S$ 的方法一样吗？

准确来说是不太一样的。

求 $g$ 的核心在筛上，可以理解为我们在用 DP 去维护埃氏筛的过程；尤其需要注意的是里面去重的那个步骤，容易忘记具体内容。

而求 $S$ 的核心其实是搜。 $S$ 的计算可以看作是搜索过程的优化，因为平时写暴力的时候也可以使用类似的方法来按照质因子搜索 $[1,n]$ 内的每个数。

关于 $S$ 的细节

既然 $S$ 本质上是搜索，那么我们就应该关注其中的剪枝小技巧。

细节一： $S$ 的边界

$S$ 搜索的范围是 $[2,n]$ ，因此任何一个搜索状态下我们都不会计算 1 的贡献。这是一个有效的剪枝，同时也要求我们必须单独计算所有质数位置的值，而不能通过某些小技巧达成这一点。

细节二：搜索下一个因子

平时搜索的一种写法是，遍历所有质数，并且允许某些质数的指数为 0。这是一个比较低效的方法，因为这个情况下我们会经过无用状态（也就是我们搜索了这个状态，但并没有产生有效贡献）。因此 $S$ 的搜索过程中就直接枚举下一个因子了。

细节三：指数剪枝

为什么 $f(p^{e+1}_i)$ 可以和 $S(\lfloor\frac n{p^e_i}\rfloor,i+1)$ 如此融洽地共存呢？为什么指数的限制为 $p_i^{e+1}\le n$ 呢？

我们可能会认为， $f(p_i^{e+1})$ 和 $S(\lfloor\frac n{p^e_i}\rfloor,i+1)$ 应该在两个和式内分别计算。但实际上，如果 $\lfloor\frac n{p_i^e}\rfloor\le p_i$ ，那么 $S(\lfloor\frac n{p_i^e}\rfloor,i+1)$ 必然为 0，我们没有必要多经过一个状态浪费时间。因此这也是一个小小的剪枝技巧。

posted @ 2020-05-14 14:06 crashed 阅读(1524) 评论(0) 编辑收藏举报

刷新页面返回顶部

登录后才能查看或发表评论，立即登录或者逛逛博客园首页

阅读排行：
· 阿里最新开源QwQ-32B，效果媲美deepseek-r1满血版，部署成本又又又降低了！
· 单线程的Redis速度为什么快？
· SQL Server 2025 AI相关能力初探
· AI编程工具终极对决：字节Trae VS Cursor，谁才是开发者新宠？
· 展开说说关于C#中ORM框架的用法！

公告

昵称： crashed
园龄： 4年11个月
粉丝： 76
关注： 18

+加关注

2025年3月

日

一

二

三

四

五

六

crashed