杜教筛学习笔记 杜教筛学习笔记
接着 Dirichlet 卷积,继续学习杜教筛。
本文中一些函数的定义见此博文
应用
通过杜教筛,我们可以快速求出某一数论函数 (f) 的前缀和,即,我们可以在低于线性的时间复杂度内求出 (S(n)=sumlimits_{i=1}^n f(i))
方法
杜教筛主要运用一个公式,通过这个公式我们建立了 (S(n)) 与 (S(frac{n}{i})) 的关系式,已知两个数论函数 (f,g),(S(n)=sumlimits_{i=1}^n f(i)),则有:
我们给个证明:
证明的关键就在于中间对于 (d,i) 枚举顺序的变化,这个证明值得细品。
那么我们有了这个公式有什么用呢?我们再来导一下:
我们可以轻松地求出 (S(n)) 的值。而 (S(lfloorfrac{n}{i} floor)) 可以通过数论分块来解决。由此看来,杜教筛的时间复杂度为低于线性复杂度。
例子
看了理论,可能还有点晕,我们来看几个例子。
求莫比乌斯函数前缀和
令 (S(n)=sumlimits_{i=1}^n mu(i) (n<2^{31}))
我们根据公式来求,我们令 (f=mu,g=1),由 Dirichlet 卷积知识,我们知道 (muast 1=varepsilon)。那么:
根据数论分块的知识,(lfloorfrac{n}{i} floor) 的取值有 (O(sqrt{n})) 个。那么时间复杂度为 (O(sumlimits_{i=1}2^{frac{1}{2^i}}) approx O(n^{frac{3}{4}}))
这明显会超时,我们需要进一步优化。其实可以这样想,一部分前缀和我们先预处理出来,另一部分再通过杜教筛来算。我们设预处理前 (k) 个,后 (n-k) 个用杜教筛做,那么此时时间复杂度是 (O(k+(n-k)^{frac{3}{4}})),由基本不等式得,(k=(n-k)^{frac{3}{4}}) 原式取最小值。我们算一下发现 (k=n^{frac{2}{3}}) 是不错的选择,此时复杂度大约为 (O(n^{frac{2}{3}}))。
//Don't act like a loser.
//You can only use the code for studying or finding mistakes
//Or,you'll be punished by Sakyamuni!!!
#include<bits/stdc++.h>
#define int long long
using namespace std;
int read() {
char ch=getchar();
int f=1,x=0;
while(ch<'0'||ch>'9') {
if(ch=='-')
f=-1;
ch=getchar();
}
while(ch>='0'&&ch<='9') {
x=x*10+ch-'0';
ch=getchar();
}
return f*x;
}
const int maxn=4e7+10;
int n,mu[maxn],p[maxn/10],cnt,s[maxn];
bool is[maxn];
map<long,long> mp;//记忆化搜索
void sieve_mu(int n) {
fill(is,is+n+1,1);
mu[1]=1;
for(int i=2;i<=n;i++) {
if(is[i]) {
p[++cnt]=i;
mu[i]=-1;
}
for(int j=1;j<=cnt;j++) {
if(p[j]*i>n) {
break;
}
is[p[j]*i]=0;
if(i%p[j]==0) {
mu[i*p[j]]=0;
break;
}
mu[i*p[j]]=-mu[i];
}
}
for(int i=1;i<=n;i++){
s[i]=s[i-1]+mu[i];
}
}
int pre(int n) {
if(n<=4e7) {
return s[n];
}
if(mp[n]) {
return mp[n];
}
int j=2,ret=1;
for(int i=2;i<=n;i=j+1) {
j=(n/(n/i));//数论分块
ret-=(j-i+1)*pre(n/i);
}
return mp[n]=ret;
}
signed main() {
sieve_mu(4e7);
int t=read();
while(t--) {
printf("%lld
",pre(read()));
}
return 0;
}
求欧拉函数前缀和
令 (S(n)=sumlimits_{i=1}^n varphi(i) (n<2^{31}))
我们根据公式来求,我们令 (f=varphi,g=1),由 Dirichlet 卷积知识,我们知道 (varphiast 1=id)。那么:
和之前一样,我们还是采取一部分直接算,一部分杜教筛的想法,时间复杂度为 (O(n^{frac{2}{3}}))
//Don't act like a loser.
//You can only use the code for studying or finding mistakes
//Or,you'll be punished by Sakyamuni!!!>
#define int long long
using namespace std;
int read() {
char ch=getchar();
int f=1,x=0;
while(ch<'0'||ch>'9') {
if(ch=='-')
f=-1;
ch=getchar();
}
while(ch>='0'&&ch<='9') {
x=x*10+ch-'0';
ch=getchar();
}
return f*x;
}
const int maxn=1e7+10;
int n,phi[maxn],p[maxn/10],cnt,s[maxn];
bool is[maxn];
map<long,long> mp;
void sieve_euler(int n) {
fill(is,is+n+1,1);
phi[1]=1;
for(int i=2;i<=n;i++) {
if(is[i]) {
p[++cnt]=i;
phi[i]=i-1;
}
for(int j=1;j<=cnt;j++) {
if(p[j]*i>n) {
break;
}
is[p[j]*i]=0;
if(i%p[j]==0) {
phi[i*p[j]]=p[j]*phi[i];
break;
}
phi[i*p[j]]=(p[j]-1)*phi[i];
}
}
for(int i=1;i<=n;i++){
s[i]=s[i-1]+phi[i];
}
}
int pre(int n) {
if(n<=1e7) {
return s[n];
}
if(mp[n]) {
return mp[n];
}
int j=2,ret=(n*(n+1)/2);
for(int i=2;i<=n;i=j+1) {
j=(n/(n/i));
ret-=(j-i+1)*pre(n/i);
}
return mp[n]=ret;
}
signed main() {
sieve_euler(1e7);
int t=read();
while(t--) {
printf("%lld
",pre(read()));
}
return 0;
}
这两个代码结合起来就可以通过洛谷上的模板题。
求莫比乌斯函数平方前缀和
令 (S(n)=sumlimits_{i=1}^n mu^2(i) (n<2^{31}))
我们继续根据公式来求,我们令 (f=mu^2,g=1),由 Dirichlet 卷积知识,我们知道 (mu^2*1=mu)。那么:
所以我们在求莫比乌斯函数平方和时,要求出莫比乌斯函数的前缀和,也是杜教筛,总时间复杂度 (O(n^{frac{2}{3}}))。
求莫比乌斯函数值平方前缀和
令 (S(n)=sumlimits_{i=1}^n (mu(i))^2 (n<2^{31}))
这个题和之前的就有所不同了。我们令 (f(i)=mu(i)^2),但是要找一个方便计算的 (g) 才行。
我们选取函数 (g(x)=[x=k^2 kin N^+])。我们来算一下 (g) 和 (f) 的狄利克雷卷积。我们会发现 (fast g=1)?!简要证明一下:
((fast g)(n)=sumlimits_{dmid n} g(d)cdot f(frac{n}{d}))。首先分类讨论,如果 (d) 是一个完全平方数,只有当 (d) 是 (n) 的所有约数中最大的完全平方数时,乘积为 (1),否则 (mu(frac{n}{d})=0)。如果 (d) 不是完全平方数,那么 (g(d)=0)。得证。
那我们的计算就简化了不少。
我们甚至不用数论分块就可以解决此题。采用同样的优化方法,时间复杂度还是约为 (O(n^{frac{2}{3}}))
相关练习题: