min_25筛入门

目录

• 1.什么是min_25筛
• 2.前置知识
  • 2.1.数论函数
  • 2.2.埃拉托色尼筛
  • 2.3.欧拉筛
• 3.min_25筛
  • 3.1.计算质数贡献
  • 3.2.计算总贡献
  • 3.3.实现
• 4.例题
  • 4.1.[LOJ]区间素数个数
  • 4.2.[LG P4213]杜教筛

1.什么是min_25筛

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

2.前置知识

2.1.数论函数

数论函数：定义域为(mathbb{N_+})子集的函数。

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

完全积性函数：若数论函数(f(x))满足(forall a,bin 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),)题目定义的(,......)

常见的完全积性函数有：

(id(x)=x; epsilon(x)=[x=1]; I(x)=1; ......)

一些有趣的性质：

[1. mu * I = epsilon ]

[2. mu * id = varphi ]

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

2.2.埃拉托色尼筛

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

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

2.3.欧拉筛

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

3.min_25筛

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

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

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

3.1.计算质数贡献

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

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

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

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

考虑这样一个 DP ：

(g(a,b))：前(a)个数进行(b)轮埃筛之后的贡献和。

可以表达为：

[g(a,b)=sum_{i=2}^a [iin P ext{或者} mp(i)> p(b)]i^k ]

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

初始状态为：

[g(a,0)=sum_{i=2}^a i^k ]

可以得到转移为：

[g(a,b)=egin{cases}g(a,b-1)& a< p(b)^2\g(a,b-1)-p(b)^k(g(lfloorfrac a{p(b)} floor,b-1) - g(p(b-1),b-1))&age p(b)^2end{cases} ]

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

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

3.2.计算总贡献

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

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

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

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

[S(a,b)=g(a,|P'|)-sum_{i=1}^{b-1}f(p(i))+sum_{i=b}^{|P'|}sum_{1le e, p(i)^{e+1}le a}left(S(lfloorfrac a{p(i)^e} floor,i+1)+f(p(i)^{e+1}) ight) ]

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

3.3.实现

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

根据整除分块理论，我们在代入(S(n,1))计算的时候，实际上(a)的取值只有(O(sqrt n))个。因此我们可以预处理这(O(sqrt n))个取值的离散化后下标。设其中一种取值为(x)，那么当(xle sqrt n)的时候，

我们直接将下标存下来；否则，由于(lfloorfrac n x floorle sqrt n)，我们将下标存在(lfloorfrac n x floor)里面。

4.例题

4.1.[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( '
' );
	return 0;
}

4.2.[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( '
' );
	}
	return 0;
}