博客两年没更了，其实博主有在打 ACM，但是因为平时用 Latex 而不是 Markdown 所以就算写了点东西也懒得搬，最近准备搬点知识点总结过来
用途
求积性函数 \(f\) 的前缀和 \(F(n)=\sum_n^{i=1} f(i)\).
复杂度 \(O(\frac{n^{\frac{3}{4}}}{\log{n}})\) (或者 \(O(n^{1-\epsilon})\), 存在 \(O(n^{\frac{2}{3}})\) 的树状数组优化版本, 但常数巨大)
使用条件
- \(f(p)\) 是低次多项式(感觉是完全积性函数的线性组合就行?)
- \(f(p^e)\) 可以快速求出
原理
Part 1
对于完全积性函数 \(f_k(i)=i^k\), 求解 \(\sum_{p_i \le n} f_k(p_i)\), 其中 \(p_i\) 是第 \(i\) 小的质数.
通过埃氏筛的方式只留下合法位置, 这样质数只需要筛到 \(\sqrt{n}\) 就可以把合数位置全部去除.
设 \(g_k(n,i)\) 为用前 \(i\) 个质数做埃氏筛, 前 \(n\) 个数剩余合法位置的 \(f_k\) 的和, 即
\[ g_k(n,i)=\sum_{j=2}^{n}{[MinP(j) > p_i \lor j \in P] \cdot f_k(j)} \]其中, \(MinP(j)\) 为 \(j\) 的最小质因数.
考虑埃氏筛的过程:
首先, \(1\) 从一开始就被筛去, 因此 \(g_k(n,0) = \sum_{i=2}^{n} f_k(p_i)\), 即自然数幂和减 \(1\).
接下来每个数被其最小质因数筛去.
\[ g_k(n,i)=g_k(n,i-1)-\sum_{j=p_i+1}^{n}{[MinP(j)=P]\cdot f_k(j)}\\ =g_k(n,i-1)-f_k(p_i)\cdot\sum_{j=p_i+1}^{n}{[MinP(j)=P]\cdot f_k(\frac{j}{p_i})}\\ =g_k(n,i-1)-f_k(p_i)\cdot\sum_{j=p_i}^{\lfloor \frac{n}{p_i} \rfloor}{[MinP(j)\ge P]\cdot f_k(j)}\\ g_k(n,i-1)-f_k(p_i)\cdot(g_k(\lfloor \frac{n}{p_i} \rfloor,i-1)-\sum_{j=1}^{i-1}f_k(p_i)) \]记 \(\sum_{j=1}^{i-1}f_k(p_i)\) 为 \(sp_k(i-1)\), 可以线性筛预处理.
根据这个式子可以dp, 但时间和空间复杂度都不能接受.
可以注意到,第 \(i\) 次筛后的 \(g\) 数组是由第 \(i-1\) 次筛后的 \(g\) 数组转移而来, 且转移点 \(\lfloor \frac{n}{p_i} \rfloor < n\), 所以可以直接省掉一维.
注意到每次转移点都是 \(\lfloor \frac{n}{x} \rfloor\) 的形式,
由于 \(\lfloor \frac{\lfloor \frac{n}{a} \rfloor}{b} \rfloor = \lfloor \frac{n}{ab} \rfloor\),
所以为了求 \(g_k(n,i)\) 实际需要用到的位置不超过 \(2\sqrt{n}\) 个.
将这 \(2\sqrt{n}\) 个位置拼成一个数组, 只计算这些位置的值, 看起来就非常优秀了.
由于对于 \(g_k(n,\infty)\), 只需要枚举不超过 \(\sqrt{n}\) 的质数, 总复杂度为:
\[ T(n)=\sum_{i=1}^{\sqrt{n}}O(\pi(\sqrt{i}))+\sum_{i=1}^{\sqrt{n}}O(\pi(\sqrt{\lfloor \frac{n}{i} \rfloor}))\\ =\sum_{i=1}^{\sqrt{n}}O(\pi(\sqrt{\frac{n}{i}}))\\ =\sum_{i=1}^{\sqrt{n}}O(\frac{\sqrt{\frac{n}{i}}}{\log(\frac{n}{i})})\\ =O(\int_{1}^{\sqrt{n}}\frac{\sqrt{\frac{n}{x}}}{\log(\sqrt{\frac{n}{x}})} \mathrm{d}x)\\ =O(\frac{\sqrt{n}}{\log n}\int_{1}^{\sqrt{n}}x^{-\frac{1}{2}} \mathrm{d}x)\\ =O(\frac{n^{\frac{3}{4}}}{\log n}) \]这样,就对所有 \(m = \lfloor \frac{n}{x} \rfloor\) 的位置都求出了\(\sum_{p_i \le m} f_k(p_i)\), 如果\(f(i)\)只在质数位置不为\(0\)(如求质数个数等, 就已经做完了).
Part 2
接下来求解 \(\sum_{i=1}^{n}f(i)\), 设:
\[S(n,i)=\sum_{j=2}^{n}{[MinP(j)>p_i] \cdot f(j)} \]所求即为 \(S(n,0)+1\).
分别考虑质数和合数位置对 \(S(n,i)\) 的贡献:
-
对质数位置, \(f(i)\) 是低次多项式, \(i^k\)的贡献为 \(g_k(n,\infty)-sp_k(i)\), 容易算出总贡献.
-
对合数位置, 需要枚举其最小质因子及其次数, 根据积性函数的性质知道,贡献为:
似乎因为某些深刻的原因(好像是只会每个位置访问一次)直接递归计算就行, 复杂度是 \(O(n^{1-\epsilon})\) 的, 但常数巨小, 在\(n \le 10^{13}\)时和\(O(\frac{n^{\frac{3}{4}}}{\log{n}})\)等数量级.
如果也写成dp递推就同样是 \(O(\frac{n^{\frac{3}{4}}}{\log{n}})\) 的.
例题
\(f(p^k)=p^k(p^k-1)\), 求 \(\sum_{i=1}^{n} f(i)\)
ll n;
int sq,m;
int tot;
int prime[N];
bool mark[N];
int sp1[N]; // f1 在前 n 个质数位置的前缀和
int sp2[N]; // f2 在前 n 个质数位置的前缀和
int g1[N*2]; // f1 在不超过 m = n/x 的质数位置的前缀和
int g2[N*2]; // f2 在不超过 m = n/x 的质数位置的前缀和
ll w[N*2];
int Ind1[N];
int Ind2[N];
// 注意 n 是全局变量
inline int& index(const ll &x){
return x<=sq?Ind1[x]:Ind2[n/x];
}
// f'(x) 与 积性函数 f(x) 在质数位置相等
// 且 f'(x) 为 f2(x) 与 f1(x) 的差
// f'(x) = x * (x - 1) = x^2 - x
inline ll f(ll x){
x%=p;
return x*(x-1)%p;
}
// f1(x) = x
inline ll f1(ll x){
return x;
}
// f2(x) = x^2
inline ll f2(ll x){
return x*x%p;
}
// f1 的前缀和
inline ll s1(ll x){
x%=p;
return x*(x+1)/2%p;
}
// f2 的前缀和
inline ll s2(ll x){
x%=p;
return x*(x+1)%p*(2*x+1)%p*inv6%p;
}
inline ll S(ll x,int j){
if(prime[j]>x)return 0;
// 质数位置的贡献
int ans=sub(sub(g2[index(x)],g1[index(x)]),sub(sp2[j],sp1[j]));
for(int i=j+1;i<=tot&&(ll)prime[i]*prime[i]<=x;++i){
for(ll e=1,pk=prime[i];pk<=x;++e,pk*=prime[i]){
// 合数位置的贡献(枚举最小质因子及其次幂)
ans=add(ans,mul(f(pk),add(S(x/pk,i),(e>1))));
}
}
return ans;
}
int main(){
r(n);
sq=sqrt(n);
for(int i=2;i<=sq;++i){
if(!mark[i]){
prime[++tot]=i;
sp1[tot]=add(sp1[tot-1],f1(i));
sp2[tot]=add(sp2[tot-1],f2(i));
}
for(int j=1,tmp;j<=tot&&(tmp=i*prime[j])<=sq;++j){
mark[tmp]=1;
if(i%prime[j]==0)break;
}
}
for(ll i=1,j;i<=n;i=j+1){
w[++m]=n/i;
j=n/w[m];
index(w[m])=m;
g1[m]=sub(s1(w[m]),1); // 先筛掉 1, f1(1)=1
g2[m]=sub(s2(w[m]),1); // 先筛掉 1, f2(1)=1
}
for(int i=1;i<=tot;++i){
for(int j=1;j<=m&&(ll)prime[i]*prime[i]<=w[j];++j){
g1[j]=sub(g1[j],mul(f1(prime[i]),sub(g1[index(w[j]/prime[i])],sp1[i-1])));
g2[j]=sub(g2[j],mul(f2(prime[i]),sub(g2[index(w[j]/prime[i])],sp2[i-1])));
}
}
printf("%d\n",add(S(n,0),1)); // 还原 f(1) 的贡献(因为是积性函数所以一定是 1 )
return 0;
}