min_25筛入门

时间:2023-11-11 12:55:32

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,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,\mu,\)题目定义的\(,......\)

        常见的完全积性函数有:

        \(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(n\log\log n)\)

        注意:对于质数\(p\),被筛掉数一定\(\ge p^2\);否则它一定会有\(< p\)质因子。

2.3.欧拉筛

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

3.min_25筛

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

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

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

3.1.计算质数贡献

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

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

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

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

        考虑这样一个 DP :

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

        可以表达为:

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

        设\(P'=[1,\lfloor\sqrt n\rfloor]\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)=\begin{cases}g(a,b-1)& a< p(b)^2\\g(a,b-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)\)。

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

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_{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)\)可以递归下去继续算。

3.3.实现

        首先需要对\([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\)里面。

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( '\n' );
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( '\n' );
}
return 0;
}