前置芝士
进入正题
现在要求一个函数 f f f 的前缀和 S S S。
但是往往连 O ( n ) O(n) O(n) 的算法都满足不了毒瘤的出题人,所以需要一个更快的做法,于是时间复杂度 O ( n 2 3 ) O(n^{\frac 2 3}) O(n32) 的杜教筛出现了。
具体做法:
找两个积性函数 g , h g,h g,h,满足 f ∗ g = h f*g=h f∗g=h。
那么有:
h
(
i
)
=
∑
d
∣
i
f
(
d
)
×
g
(
i
d
)
∑
i
=
1
n
h
(
i
)
=
∑
i
=
1
n
∑
d
∣
i
f
(
d
)
×
g
(
i
d
)
∑
i
=
1
n
h
(
i
)
=
∑
i
=
1
n
g
(
i
)
∑
j
=
1
⌊
n
i
⌋
f
(
j
)
∑
i
=
1
n
h
(
i
)
=
∑
i
=
1
n
g
(
i
)
S
(
⌊
n
i
⌋
)
h(i)=\sum_{d|i} f(d) \times g(\frac i d)\\ \sum_{i=1}^n h(i)=\sum_{i=1}^n \sum_{d|i} f(d) \times g(\frac i d)\\ \sum_{i=1}^n h(i)=\sum_{i=1}^n g(i) \sum_{j=1}^{\lfloor \frac n i \rfloor} f(j)\\ \sum_{i=1}^n h(i)=\sum_{i=1}^n g(i) S(\lfloor \frac n i \rfloor)
h(i)=d∣i∑f(d)×g(di)i=1∑nh(i)=i=1∑nd∣i∑f(d)×g(di)i=1∑nh(i)=i=1∑ng(i)j=1∑⌊in⌋f(j)i=1∑nh(i)=i=1∑ng(i)S(⌊in⌋)
把右边的第一项提出来,就有:
∑
i
=
1
n
h
(
i
)
=
S
(
n
)
g
(
1
)
+
∑
i
=
2
n
g
(
i
)
S
(
⌊
n
i
⌋
)
S
(
n
)
g
(
1
)
=
∑
i
=
1
n
h
(
i
)
−
∑
i
=
2
n
g
(
i
)
S
(
⌊
n
i
⌋
)
\sum_{i=1}^n h(i)=S(n)g(1)+\sum_{i=2}^n g(i) S(\lfloor \frac n i \rfloor)\\ S(n)g(1)=\sum_{i=1}^n h(i)-\sum_{i=2}^n g(i) S(\lfloor \frac n i \rfloor)
i=1∑nh(i)=S(n)g(1)+i=2∑ng(i)S(⌊in⌋)S(n)g(1)=i=1∑nh(i)−i=2∑ng(i)S(⌊in⌋)
那么如果 h h h 和 g g g 的前缀和很好求,那么用一个除法分块加递归就可以做了。
需要注意的是,递归里面需要对每一次求出来的前缀和记忆化,这个可以用哈希或map实现,这样时间复杂度就是 O ( n 3 4 ) O(n^{\frac 3 4}) O(n43)。
你可能会问,这样递归下去求出来的前缀和数量不会很多吗?事实上只有 2 n 2\sqrt n 2n 个,也就是不同的 ⌊ n i ⌋ \lfloor \frac n i \rfloor ⌊in⌋ 的个数,因为每一个求到的前缀和都可以表示为 ⌊ n x 1 x 2 . . . ⌋ \lfloor \frac n {x_1x_2...} \rfloor ⌊x1x2...n⌋ 的形式。
假如线性筛预处理
S
S
S 的一部分,那么上面的杜教筛的时间复杂度就可以更加优秀,比如预处理前
K
K
K 个,根据一些我不会的微积分知识可以算出复杂度会变成
O
(
K
+
n
K
)
O(K+\dfrac n {\sqrt K})
O(K+Kn),取
K
=
n
2
3
K=n^{\frac 2 3}
K=n32 时有最优复杂度,即
O
(
n
2
3
)
O(n^{\frac 2 3})
O(n32)。
实例
拿筛 μ \mu μ 和筛 φ \varphi φ 做例子。
筛莫比乌斯函数的前缀
设 f = μ f=\mu f=μ。
因为 μ \mu μ 有一个很优秀的性质—— ∑ d ∣ n μ ( d ) = [ n = 1 ] \sum_{d|n}\mu(d)=[n=1] ∑d∣nμ(d)=[n=1],写成狄利克雷卷积的形式就是: f ∗ I = ϵ f*I=\epsilon f∗I=ϵ。(不认识这些函数的话请将鼠标滚轮上滑,点开狄利克雷卷积那篇博客)
那么此时使函数 g = I , h = ϵ g=I,h=\epsilon g=I,h=ϵ,就可以做杜教筛了。
把上面的柿子拿来,就是:
S
(
n
)
g
(
1
)
=
∑
i
=
1
n
h
(
i
)
−
∑
i
=
2
n
g
(
i
)
S
(
⌊
n
i
⌋
)
S(n)g(1)=\sum_{i=1}^n h(i)-\sum_{i=2}^n g(i) S(\lfloor \frac n i \rfloor)
S(n)g(1)=i=1∑nh(i)−i=2∑ng(i)S(⌊in⌋)
显然等于:
S
(
n
)
=
1
−
∑
i
=
2
n
1
×
S
(
⌊
n
i
⌋
)
S(n)=1-\sum_{i=2}^n 1\times S(\lfloor \frac n i \rfloor)
S(n)=1−i=2∑n1×S(⌊in⌋)
代码实现:
#define maxn 6000010
int t;
ll n;
ll mu[maxn];
int prime[maxn],tt=0;
bool v[maxn];
void work()//预处理mu的前6000000位的前缀和
{
mu[1]=1ll;
for(int i=2;i<=maxn-10;i++)
{
if(!v[i])prime[++tt]=i,mu[i]=-1;
for(int j=1;j<=tt&&i*prime[j]<=maxn-10;j++)
{
v[i*prime[j]]=true;
if(i%prime[j]==0)
{
mu[i*prime[j]]=0;
break;
}
mu[i*prime[j]]=-mu[i];
}
}
for(int i=2;i<=maxn-10;i++)
mu[i]+=mu[i-1];
}
map<int,int> mapmu;
int solvemu(int x)
{
if(x<=maxn-10)return mu[x];//假如在6000000位之前
if(mapmu[x])return mapmu[x];//假如之前算过
int re=1,l=2,r;
while(l<=x)//除法分块
{
r=x/(x/l);
re-=(r-l+1)*solvemu(x/l);
l=r+1;
}
mapmu[x]=re;//记忆化
return re;
}
筛欧拉函数的前缀和
因为欧拉函数有这样一个优秀的性质: ∑ d ∣ n φ ( d ) = n \sum_{d|n}\varphi(d)=n ∑d∣nφ(d)=n,所以有: φ ∗ I = i d \varphi * I=id φ∗I=id。
于是使
g
=
I
g=I
g=I,
h
=
i
d
h=id
h=id,用上面的公式套下来就是:
S
(
n
)
=
∑
i
=
1
n
h
(
i
)
−
∑
i
=
2
n
1
×
S
(
⌊
n
i
⌋
)
S(n)=\sum_{i=1}^n h(i)-\sum_{i=2}^n 1 \times S(\lfloor \frac n i \rfloor)
S(n)=i=1∑nh(i)−i=2∑n1×S(⌊in⌋)
h h h 显然是个等差数列,可以 O ( 1 ) O(1) O(1) 求。
代码和上面类似,就不贴了。
顺手补上模板题传送门。
代码如下:
#include <cstdio>
#include <cstring>
#include <map>
#include <algorithm>
using namespace std;
#define ll long long
#define maxn 4000010
#define uint unsigned int
int t;
uint n;
ll mu[maxn],phi[maxn];
int prime[maxn],tt=0;
bool v[maxn];
void work()
{
mu[1]=phi[1]=1ll;
for(int i=2;i<=maxn-10;i++)
{
if(!v[i])prime[++tt]=i,mu[i]=-1,phi[i]=i-1;
for(int j=1;j<=tt&&i*prime[j]<=maxn-10;j++)
{
v[i*prime[j]]=true;
if(i%prime[j]==0)
{
mu[i*prime[j]]=0;
phi[i*prime[j]]=phi[i]*(ll)prime[j];
break;
}
mu[i*prime[j]]=-mu[i];
phi[i*prime[j]]=phi[i]*phi[prime[j]];
}
}
for(int i=2;i<=maxn-10;i++)
mu[i]+=mu[i-1],phi[i]+=phi[i-1];
}
map<uint,ll> mapphi;
int tot=0;
ll solvephi(uint x)
{
if(x<=maxn-10)return phi[x];
if(mapphi[x])return mapphi[x];
ll re=(ll)(1+x)*x/2ll;
uint l=2,r;
while(l<=x)
{
r=x/(x/l);
re-=(ll)(r-l+1)*solvephi(x/l);
l=r+1;
}
mapphi[x]=re;
return re;
}
map<int,ll>mapmu;
ll solvemu(uint x)
{
if(x<=maxn-10)return mu[x];
if(mapmu[x])return mapmu[x];
ll re=1;
uint l=2,r;
while(l<=x)
{
r=x/(x/l);
re-=(ll)(r-l+1)*solvemu(x/l);
l=r+1;
}
mapmu[x]=re;
return re;
}
int main()
{
work();
scanf("%d",&t);
while(t--)
{
scanf("%u",&n);
printf("%lld %lld\n",solvephi(n),solvemu(n));
}
}
但事实上这道模板题还可以做的更优。
利用
φ
\varphi
φ 和
μ
\mu
μ 之间的一个性质:
φ
=
μ
∗
i
d
\varphi = \mu * id
φ=μ∗id,那么
φ
\varphi
φ 的前缀和就可以写成:
∑
i
=
1
n
φ
(
i
)
=
∑
i
=
1
n
∑
d
∣
i
μ
(
d
)
×
i
d
=
∑
i
=
1
n
i
∑
j
=
1
⌊
n
i
⌋
μ
(
j
)
=
∑
i
=
1
n
i
×
S
(
⌊
n
i
⌋
)
\sum_{i=1}^n \varphi(i)=\sum_{i=1}^n \sum_{d|i} \mu(d) \times \frac i d\\ =\sum_{i=1}^n i \sum_{j=1}^{\lfloor \frac n i \rfloor} \mu(j)\\ =\sum_{i=1}^n i \times S(\lfloor \frac n i \rfloor)\\
i=1∑nφ(i)=i=1∑nd∣i∑μ(d)×di=i=1∑nij=1∑⌊in⌋μ(j)=i=1∑ni×S(⌊in⌋)
这样就可以用除法分块来做,转化成了求 μ \mu μ 的前缀和,因为 μ \mu μ 的前缀和我们有记忆化,所以此时跑起来飞快。
代码如下:
#include <cstdio>
#include <cstring>
#include <map>
#include <algorithm>
using namespace std;
#define maxn 5000010
#define uint unsigned int
#define ll long long
uint T,n;ll ans1,ans2;
int prime[maxn],t=0;
ll mu[maxn];
bool v[maxn];void work()
{
mu[1]=1;
for(int i=2;i<=maxn-10;i++)
{
if(!v[i])prime[++t]=i,mu[i]=-1;
for(int j=1;j<=t&&i*prime[j]<=maxn-10;j++){
v[i*prime[j]]=true;
if(i%prime[j]==0)break;
mu[i*prime[j]]=-mu[i];
}
mu[i]+=mu[i-1];
}
}
map<uint,ll>S;
ll solvemu(uint x)
{
if(x<=maxn-10)return mu[x];
if(S.count(x))return S[x];
ll re=1;
for(uint l=2,r;l<=x;l=r+1)
r=x/(x/l),re-=(r-l+1)*solvemu(x/l);
S[x]=re; return re;
}
ll solvephi(uint x)
{
ll re=0;
for(uint l=1,r;l<=x;l=r+1)
r=x/(x/l),re+=1ll*(l+r)*(r-l+1)/2*solvemu(x/l);
return re;
}
int main()
{
scanf("%u",&T);work();
while(T--){
scanf("%u",&n);
ans2=solvemu(n);
ans1=solvephi(n);
printf("%lld %lld\n",ans1,ans2);
}
}
但是,这还不是杜教筛的时间复杂度最玄学的地方。
比如说现在有一个函数
f
f
f,需要求下面这个柿子的值:
∑
i
=
1
n
f
(
⌊
n
i
⌋
)
\sum_{i=1}^n f(\lfloor \frac n i \rfloor)
i=1∑nf(⌊in⌋)
显然是一个除法分块。
但是如果 n ≤ 1 0 10 n\leq 10^{10} n≤1010, f f f 无法预处理,那么就做不了了。
考虑杜教筛。
假如写除法分块,那么里面求值的语句是长这样子的:
ans+=(r-l+1)*f(n/l);
现在的问题是 f ( n / l ) f(n/l) f(n/l) 求不了。
设 S ( n ) = ∑ i = 1 n f ( i ) S(n)=\sum_{i=1}^n f(i) S(n)=∑i=1nf(i),那么 f ( n / l ) = S ( n / l ) − S ( n / l − 1 ) f(n/l)=S(n/l)-S(n/l-1) f(n/l)=S(n/l)−S(n/l−1)。
没错,你肯定猜到了,这样转化一下之后,后面的 S S S 是个前缀和,我们用杜教筛来搞他。
于是这个 f ( n / l ) f(n/l) f(n/l)的问题就解决了。
不知道从哪里蹦出来的读者: 可是你这样的时间复杂度不是 O ( n × n 2 3 ) O(\sqrt n \times n^{\frac 2 3}) O(n×n32)的吗?更慢了!
看起来是很慢的,但是实际上,由于我们对 S ( n ) S(n) S(n) 求杜教筛后, S ( n / l ) S(n/l) S(n/l) 的前缀和也都被求出来了,于是后面用到的时候复杂度就变成了 O ( 1 ) O(1) O(1) 而不需要重新筛。
所以实际上,这样做的时间复杂度就是 O ( n 2 3 + n ) O(n^{\frac 2 3}+\sqrt n) O(n32+n) 的。
这个技巧比如这题就有用到。
另外一道题: