背景
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=1∑nF(i)=F(1)+p is a prime∑nF(p)+x is a composite∑nF(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)+p∑nF(p)+p≤n, pe≤n∑F(pe)(minp>p, i∑⌊pn⌋F(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∣)=p∑nF(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=2∑n[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,j−1) 转移到
g
(
n
,
j
)
g(n,j)
g(n,j),也就是,我们要考虑第
j
j
j 次埃氏筛筛掉了什么。
首先,经历了前面
j
−
1
j-1
j−1 次埃氏筛之后,剩下的数要么是质数,要么
m
i
n
p
>
p
j
−
1
minp >p_{j-1}
minp>pj−1。所以第
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,j−1)−pjk(g(pjn,j−1)−g(pj−1,j−1)),g(n,j−1),pj2≤npj2>n
为什么要减去
g
(
p
j
−
1
,
j
−
1
)
g(p_{j-1},j-1)
g(pj−1,j−1)?
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(pj−1,j−1)=∑i=1j−1pj,前面的
g
g
g 经历
j
−
1
j-1
j−1 次埃氏筛会把质数保留下来,而我们只需要减掉合数,所以我们要把质数去掉。
而且,
p
j
−
1
p_{j-1}
pj−1 是
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=1∑jf(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
⌊b⌊an⌋⌋=⌊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=2∑n[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+pke≤n&k>x∑f(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)=p2−p,那我们分别对它的二次项,一次项进行预处理,然后求解即可。
记得最后把
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)=p−1,最后答案再加上 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.另一个计数问题
题解戳这里喵~