杜教筛原理,实现与时间复杂度分析
引例
题目描述
给定一个正整数,求
\[ans_1=\sum_{i=1}^n\varphi(i) \]\[ans_2=\sum_{i=1}^n \mu(i) \]输入格式
本题单测试点内有多组数据。
输入的第一行为一个整数,表示数据组数 \(T\)。
接下来 \(T\) 行,每行一个整数 \(n\),表示一组询问。
输出格式
对于每组询问,输出一行两个整数,分别代表 \(ans_1\) 和 \(ans_2\)。
输入输出样例 #1
输入 #1
6 1 2 8 13 30 2333
输出 #1
1 1 2 0 22 -2 58 -3 278 -3 1655470 2
说明/提示
数据规模与约定
对于全部的测试点,保证 \(1 \leq T \leq 10\),\(1 \leq n \lt 2^{31}\)。对于全部的测试点,保证 \(1 \leq T \leq 10\),\(1 \leq n \lt 2^{31}\)。
题中要求以低于线性的时间复杂度对数论函数求和。
原理
将例写为更一般的形式:
求 $$S(n)=\sum_{i=1}^{n}f(i)$$ ,其中,\(f\) 是数论函数。
构造函数 \(g\) 并令 \(h=f*g=g*f\) ,
故,
因此,
采用数论分块加速第二个求和运算,这要求:
- 可以快速计算 \(\sum_{i=1}^nh(i)\) ;
- 可以快速计算 \(\sum_{i=l}^rg(i)\) ;
实现
实现时构造出来的 \(h,g\) 已经满足了上述要求,采用记忆化搜索计算 \(S(n)\) ,数论分块加速求和。
众所周知,\(μ∗I=ϵ\),\(φ∗I=id\) 。
-->To learn more about it
于是,带入到杜教筛公式中,
即,
考虑进一步优化,对于前 \(k\) 项使用线性筛直接计算,可以证明 \(k=O(n^{2/3})\) 最优(证明参见最后一部分)。
在实现中,一般使用 unordered_map/gp_hash_table 进行记忆化。由于常数较大,一般将预处理部分的 \(k\) 在不炸空间的情况下开得尽量大一些。
以下为C++参考实现,尽管全开了long long并且未经过卡常,常数依然中规中矩,完全可以接受。record link 。
#include <bits/stdc++.h>
#define int long long
using namespace std;
const int K=5e6;
int t,n,vh1[K],vh2[K];
bool notprime[K];
vector<int>prime;
unordered_map<int,int>s1,s2;
inline int read(){
int s=0;
bool f=0;
char c=getchar();
while(!isdigit(c)){
if(c=='-')f^=1;
c=getchar();
}
while(isdigit(c)){
s=(s<<1)+(s<<3)+(c^48);
c=getchar();
}
return f==0?s:-s;
}
void pre(){
vh1[1]=1;
vh2[1]=1;
for(int i=2;i<K;++i){
if(!notprime[i]){
prime.push_back(i);
vh1[i]=i-1;
vh2[i]=-1;
}
for(int p:prime){
if(i*p>=K)break;
notprime[i*p]=1;
vh1[i*p]=vh1[i]*vh1[p];
vh2[i*p]=-vh2[i];
if(i%p==0){
vh1[i*p]=vh1[i]*p;
vh2[i*p]=0;
break;
}
}
}
for(int i=1;i<K;++i)
vh1[i]+=vh1[i-1],vh2[i]+=vh2[i-1];
}
int S1(int n){
if(n<K)return vh1[n];
if(s1[n]!=0)return s1[n];
int res=(n+1)*n/2,l=2,r;
while(l<=n){
r=n/(n/l);
res-=S1(n/l)*(r-l+1);
l=r+1;
}
return s1[n]=res;
}
int S2(int n){
if(n<K)return vh2[n];
if(s2[n]!=0)return s2[n];
int res=1,l=2,r;
while(l<=n){
r=n/(n/l);
res-=S2(n/l)*(r-l+1);
l=r+1;
}
return s2[n]=res;
}
signed main(){
pre();
t=read();
while(t--){
n=read();
cout<<S1(n)<<" "<<S2(n)<<endl;
}
return 0;
}
时间复杂度分析
此处伪证极多,点名《算法竞赛》中的“高阶小量”伪证。
先考虑有哪些 \(S(i)\) 会被计算。
设记搜时 \(S(i)\) 会使用 \(R(i)=\{j|j=\lfloor \frac ik \rfloor,k∈\mathbb{N^+},k>1\}\) 中的元素的 \(S\) 。
下证若 \(m∈R(n)\) ,则 \(R(m)\subseteq R(n)\) 。
对于任意 \(i∈R(m)\) , \(i=\lfloor \frac mx \rfloor\) ,
又 \(m=\lfloor \frac ny \rfloor\) ,
故 \(i=\lfloor \frac {\lfloor \frac ny \rfloor}x \rfloor =\lfloor \frac n{xy} \rfloor ∈R(n)\) 。
由上知只有 \(i∈R(n)\) 才会被计算。
下令 \(T(n)\) 为计算 \(S(n)\) 的时间复杂度,并忽略计算 \(h,g\) 相关的消耗。则由数论分块相关可知 \(T(n)=\sqrt n\) 。
类似数论分块,若 \(i=\lfloor \frac nx \rfloor∈R(n)\),
则 \(i≤\sqrt n\) 或 $ i=\lfloor \frac nj \rfloor,j<\sqrt n$ 。
对于朴素的杜教筛(不使用线性筛优化),其时间复杂度为,
即为 \(O(n^{3/4})\) 。
若使用线性筛预处理前 \(k\) 项,预处理部分时间复杂度为 \(O(k)\) 。当 \(k>\sqrt n\) 时时间复杂度变为,
为求上式极值,以 \(k\) 为自变量求导,令 \(\frac {dO}{dk}=1-\frac 12 nk^{-\frac 32}=0\) ,
得到 \(k=O(n^{2/3})\) ,此时总时间复杂度为 \(O(n^{2/3})\) 。