杜教筛——省选前的学习1

BZOJ 3944 ——Sum 

题目要求给定一个数$N$，$N \leq 2^{32} - 1$，求

$$ans1 = \sum_{i = 1}^N \phi (i), ans2 = \sum_{i = 1}^N \mu (i) $$

——蛤？这个怎么做？我只知道线性筛。。。

先介绍$O(N)$的算法——线性筛。

由于Euler函数和Mobius函数都是积性函数（即对于$(a, b) = 1$, $f(ab) = f(a)f(b)$），可以利用类似贾志鹏线性筛的方法$O(n)$筛出1到$N$所有数的函数值。

#include<cstdio>
#include<cstring>
typedef long long LL;
const int N = 1000050;
LL phi[N], mobius[N];
int pr[300000];
void calc(int n) {                    //计算[1, n]的phi值及mobius值 
  memset(phi, -1, sizeof phi);
  memset(mobius, -1, sizeof mobius);
  int prnum = 0, i, j;
  phi[1] = mobius[1] = 1;
  for (i = 2; i < n; ++i) {
    if (phi[i] < 0) {
      pr[prnum++] = i;
      phi[i] = i - 1;
      mobius[i] = -1;
    }
    for (j = 0; j < prnum && (LL)i * pr[j] <= n; ++j) {
      if (i % pr[j]) {
        phi[i * pr[j]] = phi[i] * (pr[j] - 1);
        mobius[i * pr[j]] = -mobius[i];
      } else {
        phi[i * pr[j]] = phi[i] * pr[j];
        mobius[i * pr[j]] = 0;
        break;
      }
    }
  }
}
int main() {                          //O(n) - O(1)
  calc(1000000);
  for (int i = 2; i <= 1000000; ++i) {
    phi[i] += phi[i - 1];
    mobius[i] += mobius[i - 1];
  }
  long long ans1 = 0, ans2 = 0;
  int T, n;
  scanf("%d", &T);
  while (T--) {
    scanf("%d", &n);
    printf("%lld %lld", phi[n], mobius[n]);
  }
  return 0;
} 

更一般的，如果我们有任意积性函数$f(i)$，那么我们可以利用贾志鹏线性筛求出每个数的最小质因数及其幂次，从而（在p为素数时$f(p^k)$可以快速计算的情况下）线性递推出所有[1,N]的函数值。

那么，下面我们来介绍重点——杜教筛。

首先介绍狄利克雷卷积：若f，g均为$\mathbb Z \rightarrow \mathbb C$的函数（称为数论函数），定义其狄利克雷卷积为

$$(f * g)(n) = \sum_{d \mid n} f(d)g\left(\frac n d\right)$$

数论函数集与狄利克雷卷积形成群，这意味着狄利克雷卷积满足（设数论函数集为$\mathbf I$）：

1. 结合律： $\forall f, g, k \in \mathbf I, (f * g) * k = f * (g * k)$

2. 单位元： $\exists \epsilon \in \mathbf I, \forall f \in \mathbf I, f * \epsilon = \epsilon * f = f$

3. 逆元：    $\forall f \in \mathbf I, \exists f \in \mathbf I, f * g = \epsilon$

4. 封闭性： $\forall f, g \in \mathbf I, f * g \in \mathbf I$

另外，狄利克雷卷积还满足交换律：

5. 交换律： $\forall f, g \in \mathbf I, f * g = g * f$

性质2的单位函数$\epsilon (n) = \left[n = 1\right] = \begin{cases}\ 1 & n = 1\\ \;0 & n > 1\\ \end{cases}$

有了狄利克雷卷积，我们来看一看杜教筛：

假设我们要计算$S(n) = \sum_{i = 1}^N f(i)$，其中$f \in \mathbf I$。

那么，假设$g \in \mathbf I, g(1) \not= 0$， 那么

$$\begin{aligned} \sum_{i = 1}^N (f * g)(i) & = \sum_{i = 1}^N \sum_{d \mid i} g(d) f\left(\frac i d\right) \\
& = \sum_{d = 1}^N g(d) \sum_{1 \leq i \leq N, d | i} f\left(\frac i d\right) \\
& = \sum_{d = 1}^N g(d) \sum_{1 \leq i \leq \lfloor \frac n d \rfloor} f(i) \\
& = \sum_{d = 1}^N g(d) S\left(\left\lfloor \frac n d \right\rfloor\right) \end{aligned}$$

即

$$\sum_{i = 1}^N (f * g)(i)  = \sum_{d = 1}^N g(d) S\left(\left\lfloor \frac n d \right\rfloor\right)$$

$$\therefore \begin{equation} \label{recursion} g(1)S(n) = \sum_{i = 1}^N (f * g)(i)  - \sum_{i = 2}^N g(i) S\left(\left\lfloor \frac n i \right\rfloor\right)\end{equation}$$

如果我们能够快速地对$g$和$f * g$求和，那么由于所有的$\left\lfloor \frac n i \right\rfloor$只有$O(\sqrt n)$种，则对于所有$i \leq \sqrt n$ 和$i > \sqrt n$， 计算$S\left(\lfloor \frac n i \rfloor\right)$的时间复杂度为

$$O\left(\sum_{i = 1}^{\left\lfloor \sqrt n \right\rfloor} \sqrt i \right) + O\left(\sum_{i = 1}^{\left\lfloor \sqrt n \right\rfloor} \sqrt{\left\lfloor \frac n i\right\rfloor}  \right)$$

而第二项明显大一些， 它等于

$$O\left(\sum_{i = 1}^{\lfloor \sqrt n \rfloor} \lfloor \sqrt{\frac n i} \rfloor \right) \approx O\left(\int_0^{\sqrt n}\sqrt{\frac n x} dx\right) = O(n^{\frac34})$$

如果$f$为积性函数，还可以通过线性筛预处理$S(i)$的前$n^{\frac23}$项就能将复杂度降到$O\left(n^{\frac23}\right)$。

终于写完了，累死我了。蛤？还有？

考虑原题，定义$\mathbf 1: \mathbb Z \rightarrow \mathbb Z, \mathbf 1(n) = 1$，易证$\mathbf 1 * \mu = \epsilon$（貌似这就是Mobius函数的定义之一）。

代入式$(\ref{recursion})$， 得

$$ans2(n) = 1  - \sum_{i = 2}^n ans2\left(\left\lfloor \frac n i \right\rfloor\right)$$

所以只要预处理出$S$的前$n^{\frac23}$项就可以$O\left(n^{\frac23}\right)$计算了。

同理，将$\mathbf 1$和$ans1$代入$(\ref{recursion})$，得

$$ans1(n) = \frac {n(n+1)}2  - \sum_{i = 2}^n ans1\left(\left\lfloor \frac n i \right\rfloor\right)$$

 AC代码：

#include<cstdio>
#include<cstring>
#include<map>
typedef long long LL;
typedef std::map<int, LL> Map;
Map _phi, _mu;
int n;
const int N = 5000000;
LL phi[N], mu[N];
int pr[1000000];
void init(int n) {                    //计算[1, n]的phi值及mu值 
  memset(phi, -1, sizeof phi);
  memset(mu, -1, sizeof mu);
  int prnum = 0, i, j;
  phi[1] = mu[1] = 1;
  phi[0] = mu[0] = 0;
  for (i = 2; i <= n; ++i) {
    if (phi[i] < 0) {
      pr[prnum++] = i;
      phi[i] = i - 1;
      mu[i] = -1;
    }
    for (j = 0; j < prnum && (LL)i * pr[j] <= n; ++j) {
      if (i % pr[j]) {
        phi[i * pr[j]] = phi[i] * (pr[j] - 1);
        mu[i * pr[j]] = -mu[i];
      } else {
        phi[i * pr[j]] = phi[i] * pr[j];
        mu[i * pr[j]] = 0;
        break;
      }
    }
  }
  //for (i = 1; i <= 30; ++i) {
  //  printf("%d %d\n", i, phi[i]);
  //}
  for (i = 2; i <= n; ++i) {
    phi[i] += phi[i - 1];
    mu[i] += mu[i - 1];
  }
}
LL CalcPhi(LL n) {
  Map::iterator it; 
  if (n < N)
    return phi[n];
  if ((it = _phi.find(n)) != _phi.end())
    return it->second;
  LL i, last, ans = (LL)n * (n + 1) >> 1;
  for (i = 2; i <= n; i = last + 1) {
    last = n / (n / i);
    ans -= (last - i + 1) * CalcPhi(n / i);
  }
  return _phi[n] = ans;
}
LL CalcMu(LL n) {
  Map::iterator it; 
  if (n < N)
    return mu[n];
  if ((it = _mu.find(n)) != _mu.end())
    return it->second;
  LL i, last, ans = 1;
  for (i = 2; i <= n; i = last + 1) {
    last = n / (n / i);
    ans -= (last - i + 1) * CalcMu(n / i);
  }
  return _mu[n] = ans;
}
int main() {
  //freopen("sum.in", "r", stdin);
  //freopen("sum.out", "w", stdout);
  init(N - 1);
  int T;
  scanf("%d", &T);
  while (T--) {
    scanf("%d", &n);
    printf("%lld %lld\n", CalcPhi(n), CalcMu(n));
  }
  return 0;
}