杜教筛学习笔记 | fa

杜教筛能够在低于线性复杂度的时间内求出积性函数前缀和。

有两个数论函数 $f, g$ ，令

$h(x) = \sum_{d | x}{f(d) g \left( \frac xd \right)}$

则称 $h$ 是 $f$ 和 $g$ 的 Dirichlet 卷积，即 $h = f * g$ 。

两个数论函数的 Dirichlet 卷积也是数论函数。

对于一个数论函数，要求出 $S(n) = \sum_{i = 1}^n{f(i)}$ 。

构造 $g, h$ 使得 $h = f * g$ ，则

$\begin{align} \sum_{i = 1}^n{h(i)} &= \sum_{i = 1}^n{\sum_{d | 1}{g(d) f \left(\frac id \right)}} \\ &= \sum_{d = 1}^n g(d){\sum_{d | i}^n f \left(\frac id \right)} \\ &= \sum_{d = 1}^n{g(d) \sum_{i = 1}^{\lfloor \frac nd \rfloor}{f(i)}} \\ &= \sum_{d = 1}^n{g(d) S \left(\left\lfloor \frac nd \right\rfloor \right)} \\ &= g(1)S(n) + \sum_{d = 2}^n{g(d) S \left(\left\lfloor \frac nd \right\rfloor \right)} \\ S(n) &= \frac{\sum_{i = 1}^n{h(i)} - \sum_{d = 2}^n{g(d) S \left(\left\lfloor \frac nd \right\rfloor \right)}}{g(1)} \end{align}$

上面的式子实际上是给出了一个关于 $S$ 的递推式。直接这样搞是 $O(n^\frac 34)$ 的。可以线性筛出前 $n^\frac 23$ 个答案，这样总复杂度是 $O(n^\frac 23)$ 的~~，都不会证~~。

$\mathrm{Id} = \varphi * \mathbf 1 \\ S(n) = \frac{n(n + 1)}{2} - \sum_{d = 2}^n{S \left( \left\lfloor \frac nd \right\rfloor \right)} \\ \epsilon = \mu * \mathbf 1 \\ S(n) = 1 - \sum_{d = 2}^n{S \left( \left\lfloor \frac nd \right\rfloor \right)} \\ \mathrm{Id}^2 = (\varphi \times \mathrm{Id}) * \mathrm{Id} \\ \varphi = \mu * \mathrm{Id}$

目前网络上的杜教筛代码大都存在两个问题：

用 std::unordered_map 的代码在判断是否存在时直接 touch 一下来判等而不是 ::count。~~违背了 STL 的美学~~（当然，由于这个值最后总是要插入的，加了个小常数不算啥事）
判断 $S_\mu$ 是否计算过时直接和 0 判等。 $S_\mu$ 的取值在某些情况下本来就是 0，每次都要重新计算一遍相当于凭空产生常数。（对理论复杂度有无影响不会算）

两个问题很好解决，对于用 std::unordered_map 的时候这两个问题可以用 ::count 一起改掉；如果用了数组优化，第二个问题需要用其他手段记录一下。

给一下洛谷模板题的代码。

如果时空要求都不严格或是像这题一样多组询问可以用 std::unordered_map 进行记忆化。

下面的代码可以 320ms（评测鸭）通过极限数据（10 组询问，2147483638 ~ 2147483647）：

#include<bitset>
#include<iostream>
#include<unordered_map>

using ll = long long;

constexpr unsigned base = 1664510;

std::bitset<1664513> v;
int p[125641], mu[1664513];
phi[1664513];
std::unordered_map<unsigned, ll> mu_, phi_;

void sieve(unsigned N); // 普通地线筛出 1~N 的 mu 和 phi 值

ll getMu(unsigned N) {
  if (N <= base) return mu[N];
  if (mu_.count(N)) return mu_[N];
  ll ans = 0;
  for (unsigned l = 2, r, L; l <= N; l = r + 1) {
    L = N / l, r = N / L;
    ans += (r - l + 1) * getMu(L);
  }
  return mu_[N] = 1ll - ans;
}

ll getPhi(unsigned N) {
  if (N <= base) return phi[N];
  if (phi_.count(N)) return phi_[N];
  ll ans = 0;
  for (unsigned l = 2, r, L; l <= N; l = r + 1) {
    L = N / l, r = N / L;
    ans += (r - l + 1) * getPhi(L);
  }
  return phi_[N] = (1ull * N * (N + 1ull) >> 1) - ans;
}

int main() {
  unsigned T, N;

  sieve(base);
  for (std::cin >> T; T; --T) {
    std::cin >> N;
    std::cout << getPhi(N) << ' ' << getMu(N) << '\n';
  }
  return 0;
}

但是，std::unordered_map 的大常数有时令人难以接受。发现对于每次询问，需要的 $S$ 值都是形如 $S \left(\frac nd \right)$ 的，因此不会超过 $\sqrt n$ 个。

可以再开一个大小为 $\sqrt n$ 的数组进行记忆化。（我的实现中用的是 $n^\frac 23$ ）

注意对于每个 $n$ ， $f \left(\frac nd \right)$ 的值是不一样的，因此每次询问需要清空重新计算。

// 单次询问更优，但多次询问在极限数据下只能通过 T = 5
// 不过谷的板题数据很水所以可以过（

#include<bitset>
#include<iostream>

using ll = long long;

constexpr unsigned base = 1664510;

std::bitset<1664513> v, vm, vp;
int MX, p[125641]， mu[1664510], mu_[1664510];
ll phi[1664510], phi_[1664510];

inline int to(int x) {
  return MX / x;
}

void sieve(unsigned N); // 同上

int getMu(unsigned N) {
  if (N <= base) return mu[N];
  if (vm[to(N)]) return mu_[to(N)];
  int ans = 0;
  for (unsigned l = 2, r, L; l <= N; l = r + 1) {
    L = N / l, r = N / L;
    ans += (r - l + 1) * getMu(L);
  }
  vm.set(to(N));
  return mu_[to(N)] = 1 - ans;
}

ll getPhi(unsigned N) {
  if (N <= base) return phi[N];
  if (vp[to(N)]) return phi_[to(N)];
  ll ans = 0;
  for (unsigned l = 2, r, L; l <= N; l = r + 1) {
    L = N / l, r = N / L;
    ans += (r - l + 1) * getPhi(L);
  }
  vp.set(to(N));
  return phi_[to(N)] = (1ull * N * (N + 1ull) >> 1) - ans;
}

int main() {
  unsigned T, N;

  sieve(base);
  for (std::cin >> T; T; --T) {
    vm.reset(), vp.reset();

    std::cin >> N;
    MX = N;
    std::cout << getPhi(N) << ' ' << getMu(N) << '\n';
  }
  return 0;
}