Min_25篩學習小記

前言

Min_25篩的時間複雜度和洲閣篩（我不會）一樣，都是 O ( n 3 4 l o g n ) O(\frac{n^{\frac{3}{4}}}{logn}) O(lognn43​​)。但Min_25篩的時空常數，代碼複雜度據說比洲閣篩要優秀，甚至可以替代掉洲閣篩，于是我這個學不懂洲閣篩的菜雞就趕緊跑來學習了一發。

我看的是zzq的部落格

目标

設 F ( x ) F(x) F(x)是一個積性函數，現在要求 ∑ i = 1 n F ( i ) \sum_{i=1}^nF(i) i=1∑n​F(i)

Min_25篩

若要求一個積性函數的字首和，我們肯定要先把 F ( p j ) F(p_j) F(pj​)的和求出來，其中 p j p_j pj​是一個素數。

首先把不大于 n \sqrt n n

​的素數都篩出來。

由于 F ( p j ) F(p_j) F(pj​)大都是素數的幂的形式，是以考慮如何求素數幂的和。

設 g ( x , j ) = ∑ i = 2 x i k [ i 的 最 小 素 因 子 大 于 p j 或 i 是 素 數 ] g(x,j)=\sum_{i=2}^xi^k[i的最小素因子大于p_j或i是素數] g(x,j)=i=2∑x​ik[i的最小素因子大于pj​或i是素數]

若 p j 2 > x p_j^2>x pj2​>x則 g ( x , j ) = g ( x , j − 1 ) g(x,j)=g(x,j-1) g(x,j)=g(x,j−1)，否則根據定義不難得到轉移 g ( x , j ) = g ( x , j − 1 ) − p j k ( g ( ⌊ x p j ⌋ , j − 1 ) − g ( p j − 1 , j − 1 ) ) g(x,j)=g(x,j-1)-p_j^k(g(\lfloor\frac{x}{p_j}\rfloor,j-1)-g(p_j-1,j-1)) g(x,j)=g(x,j−1)−pjk​(g(⌊pj​x​⌋,j−1)−g(pj​−1,j−1))

接下來我們就可以對積性函數求和了

設 S ( x , j ) = ∑ i = 2 x F ( i ) [ i 的 最 小 素 因 子 不 小 于 p j ] S(x,j)=\sum_{i=2}^xF(i)[i的最小素因子不小于p_j] S(x,j)=i=2∑x​F(i)[i的最小素因子不小于pj​]

在求 S ( x , j ) S(x,j) S(x,j)的時候，先把所有素數的 F ( x ) F(x) F(x)求出來，這個可以通過之前處理好的 g g g數組來求。

接下來從 p j p_j pj​開始枚舉 i i i的最小素因子，設為 p k p_k pk​，再枚舉一個 e e e，滿足 p k e + 1 &lt; = x p_k^{e+1}&lt;=x pke+1​<=x，不難得到轉移 S ( x , j ) + = S ( ⌊ i p k e ⌋ , k + 1 ) F ( p k e ) + F ( p k e + 1 ) S(x,j)+=S(\lfloor\frac{i}{p_k^e}\rfloor,k+1)F(p_k^e)+F(p_k^{e+1}) S(x,j)+=S(⌊pke​i​⌋,k+1)F(pke​)+F(pke+1​)

隻要每次遞歸下去求就好了。

注意到如果每次遞歸下去求 g g g和 S S S的話，有用的 x x x隻有 O ( n ) O(\sqrt n) O(n

​)種。那麼我們隻要把第一維離散一下，然後就可以在較快的時間内求出答案。

代碼

LibreOJ #6053. 簡單的函數

#include<iostream>
#include<cstdio>
#include<cstdlib>
#include<cstring>
#include<algorithm>
#include<cmath>
using namespace std;

typedef long long LL;

const int N=200005;
const int MOD=1000000007;

int tot,s[N],g[N],id1[N],id2[N],B,m,h[N];
bool not_prime[N];
LL n,w[N],p[N];

void get_prime(int n)
{
	for (int i=2;i<=n;i++)
	{
		if (!not_prime[i]) p[++tot]=i,s[tot]=(s[tot-1]+i)%MOD;
		for (int j=1;j<=tot&&i*p[j]<=n;j++)
		{
			not_prime[i*p[j]]=1;
			if (i%p[j]==0) break;
		}
	}
}

int F(LL x,int y)
{
	if (x<=1||p[y]>x) return 0;
	int k=x<=B?id1[x]:id2[n/x],ans=(g[k]-s[y-1]-h[k]+y-1)%MOD;
	if (y==1) ans+=2;
	for (int i=y;i<=tot&&(LL)p[i]*p[i]<=x;i++)
	{
		LL t1=p[i],t2=(LL)p[i]*p[i];
		for (int e=1;t2<=x;e++,t1=t2,t2=(LL)t2*p[i])
			(ans+=((LL)F(x/t1,i+1)*(p[i]^e)%MOD+(p[i]^(e+1)))%MOD)%=MOD;
	}
	return ans;
}

int main()
{
	scanf("%lld",&n);
	B=sqrt(n);
	get_prime(B);
	for (LL i=1,last;i<=n;i=last+1)
	{
		last=n/(n/i);w[++m]=n/i;
		h[m]=(w[m]-1)%MOD;
		g[m]=(w[m]&1)?(w[m]+1)/2%MOD*(w[m]%MOD)%MOD:w[m]/2%MOD*(w[m]%MOD+1)%MOD;
		g[m]--;
		if (n/i<=B) id1[n/i]=m;
		else id2[last]=m;
	}
	for (int j=1;j<=tot;j++)
		for (int i=1;i<=m&&p[j]*p[j]<=w[i];i++)
		{
			int k=w[i]/p[j]<=B?id1[w[i]/p[j]]:id2[n/(w[i]/p[j])];
			(g[i]-=(LL)p[j]*(g[k]-s[j-1])%MOD)%=MOD;
			(h[i]-=h[k]-j+1)%=MOD;
		}
	int ans=F(n,1)+1;
	printf("%d",(ans+MOD)%MOD);
	return 0;
}