听说这玩意玩爆洲阁筛???
按照惯例,放上fj_zzq的blog
Min_25筛可以解决一类积性函数求和问题,
筛质数
假设我们现在要对n以内的质数求和,
先线筛出小于
n
\sqrt n
n的所有质数,设第i个为
p
i
p_i
pi,共有
p
0
p_0
p0个质数,
p
s
ps
ps为p的前缀和,
设函数:
S
(
x
,
j
)
=
∑
i
=
2
x
i
∗
[
i
为
质
数
或
i
的
最
小
质
因
子
大
于
p
j
]
S(x,j)=\sum_{i=2}^x i*[i为质数\ 或\ i的最小质因子大于p_j]
S(x,j)=∑i=2xi∗[i为质数 或 i的最小质因子大于pj]
设
S
u
m
(
n
)
=
(
n
+
1
)
∗
n
/
2
Sum(n)=(n+1)*n/2
Sum(n)=(n+1)∗n/2
显然,对于任意的x,
S
(
x
,
0
)
=
S
u
m
(
x
)
−
1
S(x,0)=Sum(x)-1
S(x,0)=Sum(x)−1,我们要求的Ans就是
S
(
n
,
p
0
)
S(n,p_0)
S(n,p0)
考虑转移:如果
p
j
2
>
x
p_j^2>x
pj2>x,
S
(
x
,
j
)
=
S
(
x
,
j
−
1
)
S(x,j)=S(x,j-1)
S(x,j)=S(x,j−1)
否则:
S
(
x
,
j
)
=
S
(
x
,
j
−
1
)
−
p
j
∗
(
S
(
⌊
x
p
j
⌋
,
j
−
1
)
−
p
s
j
−
1
)
S(x,j)=S(x,j-1)-p_j*(S(\lfloor \frac{x}{p_j} \rfloor,j-1)-ps_{j-1})
S(x,j)=S(x,j−1)−pj∗(S(⌊pjx⌋,j−1)−psj−1)
就这样一直递推下去即可得出 S ( n , p 0 ) S(n,p_0) S(n,p0)
不难发现这样还可以筛各种奇奇怪怪的东西,但还是不能筛任意的积性函数
积性函数求和
通过上面的方法可以求出大部分积性函数
f
(
x
)
f(x)
f(x),当x为质数时的和,现在要求的是:
∑
i
=
1
n
f
(
i
)
\sum_{i=1}^nf(i)
i=1∑nf(i)
考虑把上面的式子倒过来推,定义跟前面的一样:(可以用上面的S来预处理
G
(
x
,
p
0
)
G(x,p_0)
G(x,p0))
设
f
s
x
=
∑
i
=
1
x
f
(
p
i
)
fs_x=\sum_{i=1}^xf(p_i)
fsx=∑i=1xf(pi)
G
(
x
,
j
)
=
G
(
x
,
j
+
1
)
+
∑
k
=
1
(
G
(
⌊
x
p
j
k
⌋
,
j
+
1
)
−
f
s
j
)
∗
f
(
p
j
k
)
+
∑
k
=
2
f
(
p
j
k
)
G(x,j)=G(x,j+1)+\sum_{k=1}(G(\lfloor\frac{x}{p_j^k}\rfloor,j+1)-fs_j)*f(p_j^k)+\sum_{k=2}f(p_j^k)
G(x,j)=G(x,j+1)+k=1∑(G(⌊pjkx⌋,j+1)−fsj)∗f(pjk)+k=2∑f(pjk)
这样
A
n
s
=
G
(
n
,
0
)
Ans=G(n,0)
Ans=G(n,0)
复杂度为: O ( n 3 4 log ( n ) ) O(\frac{n^{\frac{3}{4}}}{\log(\sqrt n)}) O(log(n)n43)
可以发现对于第一个S,不用真的每次递归就直接这样递归下去,显然可以用一个for循环代替一个递归(见代码),
还可以不用递归的方式,常数还OK,
对于G也同理,
刚学一定要去做这几道题,学习正确筛法姿势,(点击传送至题解)
【LibreOJ #6053】简单的函数 //版子
【JZOJ 5683】【GDSOI2018模拟4.22】Prime //数据结构优化
【51NOD 1847】奇怪的数学题 //基本姿势
【51NOD 1965】奇怪的式子 //练手
Code
附上筛质数的代码
递归版:
LL Gg(LL n,int m)
{
if(n<N&&n<(LL)pr[m+1]*pr[m+1])return prs[n];
if(!m)return Sum(n);
for(;n<(LL)pr[m]*pr[m];--m);
LL ans=0;
for(;m;--m)ans=(ans-(Gg(n/pr[m],m-1)-prs[m-1])*pr[m])%mo;
return (ans+Sum(n))%mo;
}
非递归版:
void GS(LL n)
{
d[0]=0;
for(LL i=1,nx;i<=n;i=nx+1)
{
nx=n/(n/i);
LL t=n/i;
d[++d[0]]=t;
f[d[0]]=SUM(t,mo1)-1;
if(t<=M)id[t]=d[0];
else id1[n/t]=d[0];
}
fo(j,1,pr[0])
{
LL li=(LL)pr[j]*pr[j];
for(int i=1;d[i]>=li;++i)
{
LL t1=d[i]/pr[j];
int t=(t1>M)?id1[n/t1]:id[t1];
f[i]=(f[i]-(f[t]-prs[j-1])*(LL)pr[j])%mo1;
}
}
}