杜教筛学习笔记
杜教筛
杜教筛的基本形式
对于积性函数\(g(n)\)我们希望求他的前缀和\(S_g(n)\)，如果有另一积性函数\(f(n)\)满足\(f*g=h\)，且\(fh\)的前缀和易求，那么我们可以通过\(S_f(n) S_h(n)\)快速的求出\(S_g(n)\)。
对后面的\(S_g\left( \left\lfloor \frac{n}{d} \right\rfloor \right)\)进行数论分块就可以更快速的求出\(g\)的前缀和。
时间复杂度分析：
在求解过程中我们只需要关注\(S_g\)在\(\left\{ \left\lfloor \frac{n}{k} \right\rfloor|k\subseteq \mathbb{N}^* \right\}\)处的取值，这部分一共只有\(O\left( \sqrt{n} \right)\)种取值，因此整个时间复杂度是：
通过上面的分析我们可以知道杜教筛的时间复杂度是\(O\left( n^{\frac{3}{4}} \right)\)
如果能预处理出前\(n^{\frac{2}{3}}\)项的前缀和，复杂度可以优化到\(O\left( n^{\frac{2}{3}} \right)\)
一些例子：
-
\(\mu\)的前缀和，有\(\mu*I=e\)
所以，\(S_\mu(n)=1-\sum\limits_{i=2}^nS_\mu\left( \left\lfloor \frac{n}{i} \right\rfloor\right)\) -
\(\varphi\)的前缀和，有\(\varphi*I=id\)
所以，\(S_\varphi(n)=\frac{n\cdot (n+1)}{2}-\sum\limits_{i=2}^nS_\varphi\left( \left\lfloor \frac{n}{i} \right\rfloor\right)\)
贝尔级数
我们发现如果要用杜教筛求解积性函数前缀和问题，如何构造合适的狄利克雷卷积式子是很重要的一环，贝尔级数就是一个强有力的构造工具。
对于积性函数我们重点关注其在质数的幂处的取值，积性函数\(f(n)\)在质数\(p\)处的贝尔级数定义为\(F_p(x)=\sum\limits_{i=0}^{\infty}f(p^i)x^i\)
然后数论函数的狄利克雷卷积就可以变成一般多项式的卷积：
一些常见积性函数的贝尔级数：
- \(e\rightarrow 1\)
- \(I\rightarrow \sum\limits_{i=0}^{\infty}x^i=\frac{1}{1-x}\)
- \(id\rightarrow\sum\limits_{i=0}^{\infty}p^i\cdot x^i=\frac{1}{1-p\cdot x}\)
- \(id_k\rightarrow\sum\limits_{i=0}^{\infty}p^{k\cdot i}\cdot x^i=\frac{1}{1-p^k\cdot x}\)
- \(\mu\rightarrow\sum\limits_{i=0}^{\infty}\mu\left( p^i \right) \cdot x^i=1-x\)
- \(\varphi\rightarrow 1+\sum\limits_{i=1}^{\infty}p^{i-1}\cdot (p-1)\cdot x^i=1+\frac{p-1}{p}\cdot\sum\limits_{i=1}^{\infty}p^i\cdot x^i=\frac{1-x}{1-p\cdot x}\)
一些例子：
- 求积性函数\(f(p^k)=p^k+[k>0]\cdot (-1)^k\)的前缀和。\[\begin{aligned} &先写出f的贝尔级数F_p(x)=1+\sum\limits_{i=1}^{\infty}(p^i+(-1)^i)\cdot x^i=\frac{1+p\cdot x^2}{(1-p\cdot x)\cdot (1+x)}\\ &F_p(x)*(1-p\cdot x)\cdot (1+x)=1+p\cdot x^2,其中(1-p\cdot x)对应\mu\cdot id(点乘)(1+x)对应\mu\cdot \mu\\ &1+p\cdot x^2对应的函数满足h(n)=[\exists d,n=d^2 \wedge\mu^2(\sqrt{n})=1]\cdot \sqrt{n}\\ &令g(n)=(\mu\cdot \mu)*(\mu\cdot id),则S_f(n)=S_h(n)-\sum\limits_{d=2}^ng(d)\cdot S_f\left( \left\lfloor \frac{n}{d} \right\rfloor \right)\\ &\sum\limits_{i=1}^nh(i)=\sum\limits_{i=1}^{\sqrt{n}}i\cdot \mu^2(i)\\ &现在需要求g(n)的前缀和，注意到(1-p\cdot x)\cdot (1+x)\cdot \frac{1}{1-px}=1+x\\ &S_g(n)=S_{\mu\cdot \mu}(n)-\sum\limits_{d=2}^n \frac{d\cdot (d+1)}{2}\cdot S_g(\left\lfloor \frac{n}{d} \right\rfloor)\\ &然后需要求\mu\cdot \mu 的前缀和\\ &(\mu\cdot \mu)(n)=[n的最大平方因子=1]=\sum\limits_{d|n的最大平方因子}\mu(d)=\sum\limits_{d^2|n}\mu(d)\\ &S_{\mu\cdot \mu}(n)=\sum\limits_{i=1}^n\sum\limits_{d^2|i}\mu(d)=\sum\limits_{d=1}^{\sqrt{n}}\mu(d)\cdot \left\lfloor \frac{n}{d^2} \right\rfloor \end{aligned} \]
Powerful Number
PN的定义是不含有次数为\(1\)的质因子的正整数。即考虑正整数\(n\)的标准分解形式\(n=\prod\limits_{i=1}^{k} p_i^{a_i} ,\forall i{\kern 3pt}a_i\ge 2\)
这种形式的正整数一定可以写成\(a^2\cdot b^3\)，所以我们可以大概估计一下PN的数量，在\(\int_1^{\sqrt{n}}\sqrt[3]{\frac{n}{x^2}}dx=O\left( \sqrt{n} \right)\)级别。这样的性质会给求一些积性函数的前缀和带来便利。
假设我们要求积性函数\(f(n)\)的前缀和，我们可以先构造一个新的积性函数\(g(n)\)满足\(\forall p \subseteq prime,f(p)=g(p)\)。
此时必然存在一个积性函数\(h(n)\)满足\(f=g*h\)对于\(\forall p \subseteq prime,f(p)=g(1)\cdot h(p)+g(p)\cdot h(1)=h(p)+g(p)\)
而根据\(g\)的构造我们可以得出\(\forall p \subseteq prime,h(p)=0\)也就是说\(h\)只在PN处对答案有贡献，而PN的数量只有\(O\left( \sqrt{n} \right)\)!
接下来只需要求解\(h\)在PN处的取值，\(S_g\)可以用杜教筛计算。
考虑积性函数的取值只需要考虑其在质数的幂次处的取值：
由于是PN这里只需要考虑小于\(\sqrt{n}\)的质数，这部分暴力计算即可。
例子：定义积性函数\(f(p^k)=p^k\cdot (p^k-1)\quad p \subseteq prime\)求\(S_f(n)\)。
构造\(g(n)=n\cdot \varphi(n)\)有\(g(p)=p\cdot (p-1)=f(p)\quad p \subseteq prime\)
写出\(g\)的贝尔级数\(1+\sum\limits_{i=1}^{\infty}p^i\cdot p^{i-1}\cdot (p-1)\cdot x^i=\frac{1-p\cdot x}{1-p^2\cdot x}\)
注意到\(g*id=id_2\),所以\(S_g(n)=\frac{n\cdot (n+1)\cdot (2\cdot n+1)}{6}-\sum\limits_{d=2}^n \frac{n\cdot (n+1)}{2}\cdot S_g(\left\lfloor \frac{n}{d} \right\rfloor)\)
对于\(h\)的PN部分暴力计算。
code:
#include <cstdio>
#include <cstring>
#include <iostream>
#include <vector>
#define int long long
using namespace std;
const int mo = 1e9 + 7;
const int inv2 = (mo + 1) >> 1;
const int inv3 = (mo + 1) / 3;
const int inv6 = 1ll * inv2 * inv3 % mo;
int prime[200005], cnt = 0, g[2000005], n, sum[100005], ans = 0;
bool a[2000005];
vector<int> h[10005];
int mi(int o, int p) {
int yu = 1;
while (p) {
if (p & 1) yu = 1ll * yu * o % mo;
o = 1ll * o * o % mo;
p >>= 1;
}
return yu;
}
int S(int N) { return 1ll * N % mo * (N % mo + 1) % mo * inv2 % mo; }
int cal(int N) {
if (N <= 2000000) return g[N];
if (sum[n / N] != -1) return sum[n / N];
int yu = 1ll * N % mo * (N % mo + 1) % mo * ((N + N) % mo + 1) % mo;
yu = 1ll * yu * inv6 % mo;
for (int l = 2, r; l <= N; l = r + 1) {
r = N / (N / l);
yu -= 1ll * (S(r) - S(l - 1)) * cal(N / l) % mo;
yu = (yu % mo + mo) % mo;
}
return sum[n / N] = yu;
}
void dfs(int o, int p, int q) {
ans = (ans + 1ll * q * cal(n / o) % mo) % mo;
for (int tt = p; tt <= cnt; tt++) {
if (n / o < prime[tt] * prime[tt]) break;
int yu = o * prime[tt] * prime[tt];
for (int t = 2; yu <= n; yu *= prime[tt], t++) {
dfs(yu, tt + 1, 1ll * q * h[tt][t] % mo);
}
}
}
signed main() {
memset(sum, -1, sizeof(sum));
for (int i = 2; i <= 2000000; i++) {
if (!a[i]) {
prime[++cnt] = i;
g[i] = 1ll * i * (i - 1) % mo;
}
for (int j = 1; j <= cnt && i * prime[j] <= 2000000; j++) {
a[i * prime[j]] = 1;
if (i % prime[j] == 0) {
g[i * prime[j]] = 1ll * g[i] * prime[j] % mo * prime[j] % mo;
break;
}
g[i * prime[j]] = 1ll * g[i] * prime[j] % mo * (prime[j] - 1) % mo;
}
}
g[1] = 1;
for (int i = 2; i <= 2000000; i++) {
g[i] = (g[i - 1] + g[i]) % mo;
}
scanf("%lld", &n);
for (int i = 1; i <= cnt; i++) {
if (prime[i] > 100000) break;
h[i].push_back(1);
for (int j = prime[i], k = 1; j <= n; j *= prime[i], k++) {
int yu = 1ll * j % mo * (j % mo - 1) % mo;
for (int s = 0; s < k; s++) {
yu -= 1ll * h[i][s] * mi(prime[i], 2 * k - 2 * s - 1) % mo *
(prime[i] - 1) % mo;
yu = (yu % mo + mo) % mo;
}
h[i].push_back(yu);
}
}
dfs(1, 1, 1);
cout << ans << endl;
return 0;
}
