[笔记]杜教筛 & Powerful Number 筛
杜教筛
杜教筛的作用
杜教筛可以快速求出积性函数前缀和。如 \(\varphi\)，\(\mu\) 等。
什么是杜教筛
定义 \(f(x)\) 为一个积性函数，求 \(F(x) = \sum \limits_{i = 1}^{n} f(x)\)。
考虑构造函数 \(h, g\)，使得 \(h = f * g\)，即 \(h(n) = \sum \limits_{d | n} f(d) g(\dfrac{n}{d}) = \sum \limits_{xy = n} f(x)g(y)\)。
先求一下 \(h\) 的前缀和：
将原式裂项，得到 \(\sum \limits_{x \le n} g(x) F(\left \lfloor \dfrac{n}{x} \right \rfloor) = g(1)F(n) + \sum \limits_{2 \le x \le n} g(x) F(\left \lfloor \dfrac{n}{x} \right \rfloor)\)。
设 \(h\) 的前缀和为 \(H\)，则有 \(H(n) = g(1)F(n) + \sum \limits_{2 \le x \le n} g(x) F(\left \lfloor \dfrac{n}{x} \right \rfloor)\)。
移项得：\(g(1)F(n) = H(n) - \sum \limits_{2 \le x \le n} g(x) F(\left \lfloor \dfrac{n}{x} \right \rfloor)\)。
因此可以得到杜教筛的一些性质：
对于减号后面的部分，可以直接数论分块求，前提是可以快速求出 \(g\) 的前缀和 \(G\)。对于减号前面的部分，需要快速求出 \(h\) 的前缀和 \(H\)。那么 \(F\) 数组在等号前后都用到了怎么办啊？没事，递归 + 记忆化搜索完事。
因此可以得出杜教筛的一些性质：
-
待积 \(f\) 需要找到一组函数 \(g, h\)，满足 \(h = f * g\)。
-
\(h, g\) 的前缀和需要告诉求出。
根据上面的式子，也可以写出杜教筛的常用代码：
void calc(int n) {
if (n == 1) return f(1);
if (sum[n]) return sum[n];
int ans = H(n);
for (int l = 2, r; l <= n; r = l + 1) {
r = n / (n / l);
ans -= (r - l + 1) * calc(n / l);
} return sum[n] = ans;
}
其时间复杂度为 \(O(n ^ {\frac{3}{4}})\)。
杜教筛复杂度证明
设算法复杂度 \(T(n)\)，考虑到算法的复杂度实际为整除分块后递归的复杂度，以及递归后合并的复杂度。故有：
其中 \(O\left( \sum \limits_{j = 2}^{\sqrt{\frac{n}{i}}} T\left(\left \lfloor \dfrac{n}{ij} \right \rfloor \right) \right)\) 可以看作 \(O\) 意义下的高阶无穷小量。可以略去。
故有：
利用积分换求和的技巧，可以得到：
做一下这个积分，可以发现等于 \(O(n ^ {\frac{3}{4}})\)。
例题 1
求 \(\mu, \varphi\) 函数的前缀和。
套用刚才的式子，我们知道需要找到一个 \(g, h\)，使得 \(h = g * \mu\)。可以知道，\(g = I, h = \epsilon\) 时，有 \(\mu * I = \epsilon\)。而 \(I\) 的前缀和最好求不过了，就是 \(n\)。\(\epsilon\) 的前缀和就是 \(1\)。所以可以写出这样的代码：
unordered_map<int, int> mus;
int get_mu(int n) {
if (n == 1) return 1ll;
if (mus[n]) return mus[n];
int ans = 1ll;
for (int l = 2, r; l <= n; l = r + 1) {
r = n / (n / l);
ans -= (r - l + 1) * get_mu(n / l);
} return mus[n] = ans;
}
对于 \(\varphi\)，熟悉莫反应该知道一个式子：\(\varphi * I = id\)。所以选取 \(h = id, g = I\) 即可。其中 \(I\) 的前缀和前面已经说过，而 \(id\) 的前缀和就是等差数列求和公式，为 \(\dfrac{n (n + 1)}{2}\)。
unordered_map<int, int> phs;
int get_phi(int n) {
if (n == 1) return 1ll;
if (phs[n]) return phs[n];
int ans = n * (n + 1) >> 1;
for (int l = 2, r; l <= n; l = r + 1) {
r = n / (n / l);
ans -= (r - l + 1) * get_phi(n / l);
} return phs[n] = ans;
}
总代码如下：
#include <unordered_map>
#include <iostream>
#include <cstring>
#include <cstdio>
#define int long long
using namespace std;
int n;
unordered_map<int, int> mus;
int get_mu(int n) {
if (n == 1) return 1ll;
if (mus[n]) return mus[n];
int ans = 1ll;
for (int l = 2, r; l <= n; l = r + 1) {
r = n / (n / l);
ans -= (r - l + 1) * get_mu(n / l);
} return mus[n] = ans;
}
unordered_map<int, int> phs;
int get_phi(int n) {
if (n == 1) return 1ll;
if (phs[n]) return phs[n];
int ans = n * (n + 1) >> 1;
for (int l = 2, r; l <= n; l = r + 1) {
r = n / (n / l);
ans -= (r - l + 1) * get_phi(n / l);
} return phs[n] = ans;
}
signed main() {
int T; scanf("%d", &T);
while (T -- ) {
mus.clear();
phs.clear();
scanf("%lld", &n);
printf("%lld %lld\n", get_phi(n), get_mu(n));
} return 0;
}
欸，怎么只有 \(30\) 分？说明算法可以继续优化。
算法优化
假设假设我们已经用线性筛筛出了前 \(k\) 项函数前缀和，即预处理了 \(F(1) \sim F(k)\)。那么对于 \(\left \lfloor \dfrac{n}{i} \right \rfloor \le k\) 的部分都不需要递归计算了。从而新算法复杂度（这里指递归计算的部分）\(T'(n) = O\left( \sum \limits_{i = 2}^{\frac{n}{k}} \sqrt{\left \lfloor \dfrac{n}{i} \right \rfloor} \right) = O\left( \int_{i = 0}^{\frac{n}{k}} \sqrt{\left \lfloor \dfrac{n}{i} \right \rfloor} \right) = O\left( \dfrac{n}{\sqrt{k}}\right)\)
所以算法时间复杂度就是 \(T(n) = O(k) + T'(n) = O(k) + O\left( \dfrac{n}{\sqrt{k}}\right)\)
根据均值不等式可得，当 \(k = \dfrac{n}{\sqrt{k}}\)，即 \(k = n ^ {\frac{2}{3}}\) 时，时间复杂度最优为 \(O(n ^ {\frac{2}{3}})\)。
因此对于模板题，只需要筛出约前 \(6000000\) 个 \(\varphi\) 和 \(\mu\) 的前缀和就可以了。
#include <unordered_map>
#include <iostream>
#include <cstring>
#include <cstdio>
#define int long long
using namespace std;
const int N = 6000010;
int n, lim;
int p[N], mu[N], phi[N];
int smu[N], sphi[N], cnt;
bool is_prime[N];
unordered_map<int, int> mus;
void init(int n) {
mu[1] = phi[1] = 1;
for (int i = 2; i <= n; i ++ ) {
if (!is_prime[i]) p[ ++ cnt] = i, phi[i] = i - 1, mu[i] = -1;
for (int j = 1; j <= cnt and p[j] * i <= n; j ++ ) {
is_prime[i * p[j]] = true;
if (i % p[j] == 0) {
phi[i * p[j]] = phi[i] * p[j];
break;
}
mu[i * p[j]] = - mu[i];
phi[i * p[j]] = phi[i] * phi[p[j]];
}
}
for (int i = 1; i <= n; i ++ )
sphi[i] = sphi[i - 1] + phi[i],
smu[i] = smu[i - 1] + mu[i];
}
int get_mu(int n) {
if (n <= lim) return smu[n];
if (mus[n]) return mus[n];
int ans = 1ll;
for (int l = 2, r; l <= n; l = r + 1) {
r = n / (n / l);
ans -= (r - l + 1) * get_mu(n / l);
} return mus[n] = ans;
}
unordered_map<int, int> phs;
int get_phi(int n) {
if (n <= lim) return sphi[n];
if (phs[n]) return phs[n];
int ans = n * (n + 1) >> 1;
for (int l = 2, r; l <= n; l = r + 1) {
r = n / (n / l);
ans -= (r - l + 1) * get_phi(n / l);
} return phs[n] = ans;
}
signed main() {
int T; scanf("%d", &T);
lim = 6000000;
init(lim);
while (T -- ) {
mus.clear();
phs.clear();
scanf("%lld", &n);
printf("%lld %lld\n", get_phi(n), get_mu(n));
} return 0;
}
Powerful Number 筛
Powerful Number
简称 \(\text{PN}\)。\(\text{PN}\) 是指正整数 \(n\)，满足 \(n\) 的分解 \(\prod p_i ^ {c_i}\) 中，\(\forall i \in [1, k], c_i > 1\)。
暂且假设 \(\text{Powerful Number}\) 组成的集合为 \(\mathbb{NP}\)。
引理 \(1\)：\(\forall n \in \mathbb{NP}\)，\(n\) 可以表示成 \(a ^ 2 b ^ 3\) 的形式。
证明很简单，把 \(n\) 分解成 \(\prod p_i ^ {c_i}\) 形式后，如果 \(c_i\) 是偶数就归到 \(a ^ 2\) 类里，否则就减三再扔到 \(a ^ 2\) 类里面去。其本质是任意 \(\ge 2\) 的正整数都可以表示成 \(2a + 3b\) 的形式。
引理 \(2\)：\(n\) 以内的 \(\text{PN}\) 数量为 \(O(\sqrt{n})\) 级别。
证明：转化成积分形式就是 \(\int _ 0 ^ {\sqrt{n}} \left ( \dfrac{n}{x ^ 2} \right ) ^ {\frac{1}{3}} \mathrm{dx}\)。
积出来以后是 \(O(\sqrt{n})\) 级别的。
Powerful Number 筛
下面的函数均指积性函数。
对于 \(f(x)\) 求前缀和，记为 \(F(x)\)。\(F(n) = \sum \limits_{i = 1}^{n} f(i)\)。
考虑构造函数 \(g, h\)，使得 \(f = g * h\)，且使 \(g\) 在 \(x \in \mathbb{P}\) 的时候取值与 \(f(x)\) 相同。
考虑当 \(n \in \mathbb{P}\) 时，\(f(n) = \sum \limits_{d | n} g(d) h (\dfrac{n}{d}) = g(1)h(n) + g(n)h(1)\)。由于 \(g(n) = f(n), g(1) \ne 0\)，所以 \(h(n) = 0\)。由于 \(h\) 为积性函数，可以得到当 \(n \notin \mathbb{NP}\) 时，\(h(n)\) 均等于 \(0\)。
所以继续化简式子:
也就是说，我们只要能够快速算出 \(n\) 以内的所有 \(\text{Powerful Number}\)，并且能够快速计算 \(G\) 的前缀和即可。
其中 \(G\) 的前缀和可以直接杜教筛做到 \(O(n ^ {\frac{2}{3}})\) 的复杂度，而 \(h\) 的点值可以直接爆搜所有 \(\text{Powerful Number}\)。
算法复杂度为 \(O\left(\dfrac{n ^ {\frac{3}{4}}}{\log n} \right)\)
例题
\(f(p ^ k) = p ^ k (p ^ k - 1)\)
考虑设计 \(g(p) = p \times \varphi(p) = p \times (p - 1) = f(p)\)。
所以 \(g(x) = x \times \varphi(x)\)。
再设计一个 \(h(x)\)，使得 \(f = g * h\)。不想写了，反正就是 \(h(p ^ k) = (k - 1)(p - 1)p ^ k\)。
#include <unordered_map>
#include <iostream>
#include <cstring>
#include <cstdio>
#define int long long
using namespace std;
const int N = 3000010;
const int mod = 1e9 + 7;
unordered_map<int, int> Gsum;
int phi[N], p[N], s[N];
int lim, ans, n, cnt, inv2, inv6;
bool is_prime[N];
int Mod(int a, int b) {
return ((a % mod) * (b % mod)) % mod;
}
int power(int a, int b = mod - 2) {
int ans = 1;
for (; b; b >>= 1, a = a * a % mod)
if (b & 1) ans = ans * a % mod;
return ans;
}
void init(int n) {
phi[1] = 1;
for (int i = 2; i <= n; i ++ ) {
if (!is_prime[i]) p[ ++ cnt] = i, phi[i] = i - 1;
for (int j = 1; j <= cnt and i * p[j] <= n; j ++ ) {
is_prime[i * p[j]] = true;
if (i % p[j] == 0) { phi[i * p[j]] = phi[i] * p[j]; break; }
phi[i * p[j]] = phi[i] * phi[p[j]];
}
}
for (int i = 1; i <= n; i ++ )
s[i] = (s[i - 1] + i * phi[i] % mod) % mod;
}
int get_G(int n) {
if (n <= lim) return s[n];
if (Gsum[n]) return Gsum[n];
int ans = Mod(Mod(Mod(n, n + 1), (n * 2 % mod + 1)), inv6);
for (int l = 2, r; l <= n; l = r + 1) {
r = n / (n / l);
ans -= Mod(Mod(Mod(l % mod + r % mod, r - l + 1), inv2), get_G(n / l));
if (ans < 0) ans += mod;
}
return Gsum[n] = ans;
}
void dfs(int now, int num, int val) {
if (now > cnt or num * p[now] > n) {
if (n / num <= lim) ans = (ans + Mod(val, s[n / num])) % mod;
else ans = (ans + Mod(val, Gsum[n / num])) % mod;
return;
}
dfs(now + 1, num, val);
int u = Mod(Mod((p[now] - 1), p[now]), p[now]), tmp = p[now] * p[now];
for (int i = 2; num * tmp <= n; i ++ , tmp = tmp * p[now]) {
dfs(now + 1, num * tmp, Mod(Mod(val, i - 1), u));
u = Mod(u, p[now]);
}
}
signed main() {
inv2 = power(2);
inv6 = power(6);
scanf("%lld", &n);
lim = 3000000;
init(lim); get_G(n);
dfs(1, 1, 1);
printf("%lld\n", ans);
return 0;
}
