【学习笔记】min_25筛

背景

GDCPC2024
出题人:出这道 min25 筛是给大家增加过题数的 [呲牙][大哭][呲牙][大哭]

min25筛是干啥的

快速求一个积性函数 F ( x ) F(x) F(x) 的前缀和
这个 F ( x ) F(x) F(x) 需要满足:
F ( p ) = ∑ i = 0 a i p i F(p)=\sum_{i=0}a_ip^i F(p)=i=0aipi,也就说在质数处 F ( p ) F(p) F(p) 是关于 p p p 的多项式
并且 F ( p k ) F(p^k) F(pk) 需要能够快速计算
比杜教筛适用范围更广。
时间复杂度 O ( n 3 4 l o g n ) O(\frac{n^{\frac{3}{4}}}{logn}) O(lognn43)

Step 1 分类

先对 i i i 按质数和合数分类。
∑ i = 1 n F ( i ) = F ( 1 ) + ∑ p   i s   a   p r i m e n F ( p ) + ∑ x   i s   a   c o m p o s i t e n F ( x ) \sum_{i=1}^nF(i)=F(1) + \sum_{p\ is\ a\ prime}^nF(p)+\sum_{x\ is\ a\ composite}^nF(x) i=1nF(i)=F(1)+p is a primenF(p)+x is a compositenF(x)
由积性函数的性质,可以在合数那里提出一个 F ( p e ) F(p^e) F(pe)
F ( 1 ) + ∑ p n F ( p ) + ∑ p ≤ n ,   p e ≤ n F ( p e ) ( ∑ m i n p > p ,   i ⌊ n p ⌋ F ( i ) ) F(1)+\sum_{p}^nF(p)+\sum_{p\leq\sqrt{n},\ p^e\leq n}F(p^e)(\sum_{minp>p,\ i}^{\left\lfloor\frac{n}{p}\right\rfloor}F(i)) F(1)+pnF(p)+pn , penF(pe)(minp>p, ipnF(i))

Step2 预处理

g ( n , j ) g(n,j) g(n,j) 为前 n n n F ( i ) F(i) F(i) 经历 j j j 次埃氏筛,剩下的 F ( i ) F(i) F(i) 值之和。最后我们需要用到的值是
g ( n , ∣ P ∣ ) = ∑ p n F ( p ) g(n,|P|) =\sum_{p}^nF(p) g(n,P)=pnF(p)
由于 F ( p ) F(p) F(p) 是关于 p p p 的多项式,所以我们只需要令 F ( i ) = i k F(i)=i^k F(i)=ik 对应多项式中其中一项,然后再把每一项加起来。即:
g ( n , j ) = ∑ i = 2 n [ i   i s   a   p r i m e   o r   m i n p ( i ) > p j ] i k g(n,j) = \sum_{i=2}^n[i\ is\ a\ prime\ or\ minp(i)>p_j]i^k g(n,j)=i=2n[i is a prime or minp(i)>pj]ik
其中 m i n p ( i ) minp(i) minp(i) i i i 的最小质因子, p j p_j pj为第 j j j 大的质数。
那么,这玩意怎么求呢?注意,这里的 n n n 很大,不能直接线性筛。

考虑如何从 g ( n , j − 1 ) g(n,j-1) g(n,j1) 转移到 g ( n , j ) g(n,j) g(n,j),也就是,我们要考虑第 j j j 次埃氏筛筛掉了什么。
首先,经历了前面 j − 1 j-1 j1 次埃氏筛之后,剩下的数要么是质数,要么 m i n p > p j − 1 minp >p_{j-1} minp>pj1。所以第 j j j 次筛掉的数一定是 m i n p = p j minp=p_j minp=pj 的合数。
并且, p j p_j pj 能筛掉的最小的数是 p j 2 p_j^2 pj2
我们可以从这些合数提出一个 p j p_j pj,根据完全积性函数的性质,得到转移式:
g ( n , j ) = { g ( n , j − 1 ) − p j k ( g ( n p j , j − 1 ) − g ( p j − 1 , j − 1 ) ) , p j 2 ≤ n g ( n , j − 1 ) , p j 2 > n g(n,j)=\begin{cases} g(n,j-1)-p_j^k(g(\frac{n}{p_j},j-1)-g(p_{j-1},j-1)), & p_j^2 \leq n \\ g(n,j-1), & p_j^2>n \end{cases} g(n,j)={g(n,j1)pjk(g(pjn,j1)g(pj1,j1)),g(n,j1),pj2npj2>n
为什么要减去 g ( p j − 1 , j − 1 ) g(p_{j-1},j-1) g(pj1,j1)?
g ( p j − 1 , j − 1 ) = ∑ i = 1 j − 1 p j g(p_{j-1},j-1)=\sum_{i=1}^{j-1}p_j g(pj1,j1)=i=1j1pj,前面的 g g g 经历 j − 1 j-1 j1 次埃氏筛会把质数保留下来,而我们只需要减掉合数,所以我们要把质数去掉。
而且, p j − 1 p_{j-1} pj1 n \sqrt n n 以内的数,可以使用线性筛预处理前缀和,也就是
s p j = ∑ i = 1 j f ( p i ) = g ( p j , j ) sp_j=\sum_{i=1}^jf(p_i)=g(p_j,j) spj=i=1jf(pi)=g(pj,j)
那么我们要的东西 g ( n , ∣ P ∣ ) g(n,|P|) g(n,P) 的质数集合 P P P 为不超过 n \sqrt n n 的质数集。
但是 n n n 还是太大了,我们依旧无法一个一个地把 g g g 求出来。但是我们又会想到另一个重要结论:
⌊ ⌊ n a ⌋ b ⌋ = ⌊ n a b ⌋ \left\lfloor \frac{\left\lfloor \frac{n}{a}\right\rfloor}{b}\right\rfloor = \left\lfloor\frac{n}{ab}\right\rfloor ban=abn
也就说,我们需要的数都能被表示成 x x x 或者 ⌊ n x ⌋ \left\lfloor\frac{n}{x}\right\rfloor xn 的形式。这些数一共只有 O ( n ) O(\sqrt n) O(n ) 个。
那么如何访问这 O ( n ) O(\sqrt n) O(n ) 个位置呢?
如果偷懒用 m a p map map 的话会多一个 l o g log log,我们可以用两个数组 i d 1 ,   i d 2 id1,\ id2 id1, id2 存储它们的下标。

if(w[tot]<=sqr) 
	id1[w[tot]]=tot;
else 
	id2[n/w[tot]]=tot;

step3 求答案

我们令 S ( n , j ) S(n,j) S(n,j) 2 2 2~ n n n m i n p > p j minp>p_j minp>pj 的函数值之和,即:
S ( n , j ) = ∑ i = 2 n [ m i n p ( i ) > p j ] S(n,j)=\sum_{i=2}^n[minp(i)>p_j] S(n,j)=i=2n[minp(i)>pj]
那么我们要求的答案就是 S ( n , 0 ) S(n,0) S(n,0)

同样的,我们依旧把和分成质数的和,合数的和。枚举它们最小的因子,那么就有下面的转移:
S ( n , j ) = g ( n , ∣ P ∣ ) − s p j + ∑ p k e ≤ n & k > x f ( p k e ) ( S ( n p k e , k ) + [ e ≠ 1 ] ) S(n,j)=g(n,|P|)-sp_j+\sum_{p_k^e\leq n \&k>x}f(p_k^e)(S(\frac{n}{p_k^e},k)+[e≠1]) S(n,j)=g(n,P)spj+pken&k>xf(pke)(S(pken,k)+[e=1])
前面的 g ( n , ∣ P ∣ ) − s p j g(n,|P|)-sp_j g(n,P)spj 显然是质数的答案部分。
后面的合数同样可以根据积性函数的性质提出 f ( p k e ) f(p_k^e) f(pke)
需要注意的是, S ( n p k e , k ) S(\frac{n}{p_k^e},k) S(pken,k) 是会把 p k e [ e > 1 ] p_k^e[e>1] pke[e>1] 筛掉的,所以要把它给补上,也就是加上 [ e ≠ 1 ] [e≠1] [e=1]

至此,min25筛的流程就讲完啦。

例题

【模板】Min_25 筛

在这里插入图片描述
在这里插入图片描述

思路

显然, f ( p ) = p 2 − p f(p)=p^2-p f(p)=p2p,那我们分别对它的二次项,一次项进行预处理,然后求解即可。
记得最后把 f ( 1 ) f(1) f(1) 加上

代码

#include<bits/stdc++.h>
#define int long long
using namespace std;
const int N=2e5+7,inf=1e18,mod=1e9+7;
vector<int> g1(N),g2(N),id1(N),id2(N),sp1(N),sp2(N),p,w(N);
int sqr,tot=0,n;
int power(int x,int t)
{
	int b=1;
	while(t)
	{
		if(t&1) b=b*x%mod;
		x=x*x%mod; t>>=1;
	}
	return b;
}
void init(int n)
{
	vector<bool> bz(n+1);
	p.push_back(0);
	tot=0;
	bz[1]=1;
	for(int i=1; i<=n; i++)
	{
		if(!bz[i])
		{
			p.push_back(i);
			int now=p.size()-1;
			sp1[now]=(sp1[now-1]+i)%mod;
			sp2[now]=(sp2[now-1]+i*i%mod)%mod;
		}
		for(auto j:p)
		{
			if(!j) continue;
			if(i*j>n) break;
			bz[i*j]=1;
			if(i%j==0) break;
		}
	}
}
int S(int x,int y)
{
	if(x<=1||p[y]>=x) return 0;
	int k=(x<=sqr)?id1[x]:id2[n/x];
	//prime part
	int ans=(g2[k]-g1[k])-(sp2[y]-sp1[y]);
	ans%=mod;
	//other parts
	for(int k=y+1; k<p.size()&&p[k]*p[k]<=x; k++)
	{
		int t=p[k];
		for(int e=1; t<=x; e++,t*=p[k])
		{
			int p=t%mod;// avoid overflow
			(ans+=p*(p-1)%mod*(S(x/t,k)+(e!=1))%mod)%=mod;
		}
	}
	return ans;
}
void O_o()
{
	//f(p^k) = p^k(p^k - 1)
	cin>>n;
	sqr=sqrt(n);
	init(sqr);
	//let g(n,0) = \sum_{i=1}^n f'(i) 
	int inv2=power(2,mod-2),inv3=power(3,mod-2);
	for(int i=1,j; i<=n; i=j+1)
	{
		j=n/(n/i);
		w[++tot]=n/i;
		int now=w[tot]%mod;
		g1[tot]=now*(now+1)%mod*inv2%mod-1;
		g2[tot]=now*(now+1)%mod*(2*now+1)%mod*inv2%mod*inv3%mod-1;// 4+9+16...
		if(w[tot]<=sqr) id1[w[tot]]=tot;
		else id2[n/w[tot]]=tot;
	}
	//get g(n,|P|)
	for(int i=1; i<p.size(); i++)
	{
		for(int j=1; j<=tot&&p[i]*p[i]<=w[j]; j++)
		{
			//get g(w[j],i)
			int k=w[j]/p[i]<=sqr?id1[w[j]/p[i]]:id2[n/(w[j]/p[i])];
			//g(w[j],i) = g(w[j],i-1) - f'(p[i])*(g(w[k],j-1)-sp[i-1])
			(g1[j]-=p[i]*(g1[k]-sp1[i-1]+mod)%mod)%=mod;
			(g2[j]-=p[i]*p[i]%mod*(g2[k]-sp2[i-1]+mod)%mod)%=mod;
		}
	}
	int ans=S(n,0)+1;
	(ans+=mod)%=mod;
	cout<<ans<<"\n";
}
signed main()
{
	ios::sync_with_stdio(false); cin.tie(0); cout.tie(0);
	cout<<fixed<<setprecision(2);
	int T=1;
//	cin>>T;
	while(T--)
	{
		O_o();
	}
}

简单的函数

在这里插入图片描述
在这里插入图片描述

思路

先考虑质数处的值,显然除了 2 以外,所有质数都是奇数。但我们可以先令 f ( p ) = p − 1 f(p)=p-1 f(p)=p1,最后答案再加上 2。
那么我们需要对 p p p − 1 -1 1 进行预处理,然后直接套 min25 筛即可。

代码

#include<bits/stdc++.h>
#define int long long
using namespace std;
const int N=2e5+7,inf=1e18,mod=1e9+7;
vector<int> g1(N),g0(N),sp1(N),id1(N),id2(N),p,w(N);
int n,sqr,tot;
void init(int n)
{
	p.push_back(0);
	tot=0;
	vector<bool> bz(n+1);
	for(int i=2; i<=n; i++)
	{
		if(!bz[i])
		{
			p.push_back(i);
			int now=p.size()-1;
			sp1[now]=(sp1[now-1]+i)%mod;
		}
		for(auto j:p)
		{
			if(!j) continue;
			if(i*j>n) break;
			bz[i*j]=1;
			if(i%j==0) break;
		}
	}
}
int S(int x,int y)
{
	if(x<=1||p[y]>=x) return 0;
	int k=(x<=sqr)?id1[x]:id2[n/x];
	int ans=g1[k]-g0[k]-sp1[y]+y;
	ans%=mod;
	for(int k=y+1; k<p.size()&&p[k]*p[k]<=x; k++)
	{
		int t=p[k];
		for(int e=1; t<=x; e++,t*=p[k])
		{
			(ans+=(p[k]^e)*(S(x/t,k)+(e!=1))%mod)%=mod;
		}
	}
	return ans;
}
void O_o()
{
	cin>>n;
	if(n==1)
	{
		cout<<1<<"\n";
		return;
	}
	sqr=sqrt(n);
	init(sqr);
	for(int i=1,j; i<=n; i=j+1)
	{
		j=n/(n/i);
		w[++tot]=n/i;
		int now=w[tot]%mod;
		g1[tot]=now*(now+1)/2%mod-1;
		g0[tot]=now-1;
		if(w[tot]<=sqr) id1[w[tot]]=tot;
		else id2[n/w[tot]]=tot;
	}
	for(int i=1; i<p.size(); i++)
	{
		for(int j=1; j<=tot,p[i]*p[i]<=w[j]; j++)
		{
			int k=w[j]/p[i]<=sqr?id1[w[j]/p[i]]:id2[n/(w[j]/p[i])];
			//g(w[j],i) = g(w[j],i-1) - f'(p[i])*(g(w[k],j-1)-sp[i-1])
			(g1[j]-=p[i]*(g1[k]-sp1[i-1])%mod)%=mod;
			(g0[j]-=(g0[k]-(i-1))%mod)%=mod;
		}
	}
	int ans=S(n,0)+3;
	(ans+=mod)%=mod;
	cout<<ans<<"\n";
}
signed main()
{
	ios::sync_with_stdio(false); cin.tie(0); cout.tie(0);
	cout<<fixed<<setprecision(2);
	int T=1;
//	cin>>T;
	while(T--)
	{
		O_o();
	}
}

【GDCPC2024】J.另一个计数问题

题解戳这里喵~

  • 14
    点赞
  • 22
    收藏
    觉得还不错? 一键收藏
  • 0
    评论
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值