筛法的进化


筛法

埃氏(Eratosthenes,埃拉托斯特尼)筛法

考虑这样一个情景:如果 x 是合数,那么 x 的倍数也一定是合数。利用这个结论,可以避免很多次不必要的素数检测。

如果从小到大考虑每个数,然后同时把当前这个数所有比自己大的倍数标记为合数,那么没有被标记的数就是素数了。

可以参考如下程序:

 int Eratosthenes(int n)
 {
  int p = 0;
  for (int i = 0; i <= n; ++i) is_prime[i] = 1;
  is_prime[0] = is_prime[1] = 0;
  for (int i = 2; i <= n; ++i)
    {
  if (is_prime[i]) {
  prime[p++] = i;
  if ((long long)i * i <= n)
  for (int j = i * i; j <= n; j += i)
    // 因为从 2 到 i - 1 的倍数我们之前筛过了,这里直接从 i 的倍数开始,提高了运行速度
  is_prime[j] = 0;  // 是i的倍数的均不是素数
  }
  }
  return p;
 }

上面的程序就是埃氏筛法。

接下来讨论该算法的时间复杂度:

每一次对数组的操作花费一个单位时间,则时间复杂度为:O\left(n\sum^{\pi(n)}_{k = 1} {\frac{1}{p_k}}\right),其中 p_k??? 表示第 k????? 小的质数。下面进行对其进行简单的化简:

\begin{align*}\sum_{p \le x}{\frac{1}{p}} &= \frac{1}{2} + \int^x_2{\frac{\text{d} g(t)}{\log t}}\\ &= \frac{g(x)}{\log x} + \int^x_2{\frac{g(t)}{t \log^2 t}\text{d}t} \\&= 1 + \frac{h(x)}{\log x}+\int^x_2{\frac{\text{d}t}{t \log t}} + \int^\infin_2{\frac{h(t)}{t \log^2 t}\text{d}t} - \int^\infin_x{\frac{h(t)}{\log t}\text{d}t}\\ &=\log \log x + 1 - \log \log 2 - \int^\infin_2{\frac{h(t)}{t \log^2 t}\text{d}t} + O\left(\frac{1}{\log x}\right)\end{align*}?

其中 \begin{align*}1 - \log \log 2 - \int^\infin_2{\frac{h(t)}{t \log^2 t}\text{d}t}\end{align*}? 易证其收敛于 \begin{align*}B_1 = \lim_{x \to \infin} \left[\sum_{p \le x} \frac{1}{p} - \log \log x\right]\end{align*}

易得素数的倒数和可以写为: \begin{align*}\sum_{p \le x}{\frac{1}{p}} = \log \log x + B_1 + O(\frac{1}{\log x})\end{align*}?。

下面给出上述证明的简化版:

根据 \pi(n) = \Theta(\frac{n}{\log n})? 的定义,可以知道第 n 个素数的大小为: \Theta(n \log n) ,于是有:

\begin{align*}\sum^{\pi(n)}_{k = 1} {\frac{1}{p_k}} &= O\left(\sum^{\pi(n)}_{k = 2} {\frac{1}{k \log k}} \right)\\ &= O\left(\int_2^{\pi(n)} {\frac{\text{d}x}{x \log x}}\right) \\ &= O\left(\log \log \pi(n) \right) = O(\log \log n)\end{align*}??

所以埃氏筛法的时间复杂度为:\begin{align*}O\left(n \log \log n \right)\end{align*}? 。?????

但是时间复杂度仍有一些高,进行一些优化:

筛至平方根

显然,要找到至多至 n 的素数,只要找到不超过 \sqrt{n} 的素数就可以了。这种优化不会影响渐进时间复杂度,实际上重复以上证明,可以得到 n \ln \ln \sqrt{n} + o(n)?,根据对数的性质,它们的渐进相同,但操作次数会明显减少。下面给出该优化后的参考程序:

 int n;
 vector<char> is_prime(n + 1, true);
 is_prime[0] = is_prime[1] = false;
 for (int i = 2; i * i <= n; i++)
 {
  if (is_prime[i])
    {
  for (int j = i * i; j <= n; j += i)  is_prime[j] = false;
  }
 }

 

只筛奇数

因为除 2 以外的偶数都是合数,所以我们可以直接跳过它们,只用关心奇数就好。

首先,这样做能让我们内存需求减半;其次,所需的操作大约也减半。

线性筛法(欧式筛法)

埃氏筛法仍有优化空间,它会将一个合数重复多次标记。有没有什么办法省掉无意义的步骤呢?答案是肯定的。

如果能让每个合数都只被标记一次,那么时间复杂度就可以降到 O\left(n\right)??? 了。这个算法就称为“线筛”,即线性筛法。参考代码如下所示:

 void init() {
  phi[1] = 1;
  for (int i = 2; i < MAXN; ++i) {
  if (!vis[i]) {
  phi[i] = i - 1;
  pri[cnt++] = i;
  }
  for (int j = 0; j < cnt; ++j) {
  if (1ll * i * pri[j] >= MAXN) break;
  vis[i * pri[j]] = 1;
  if (i % pri[j])
  {
  phi[i * pri[j]] = phi[i] * (pri[j] - 1);
  }
  else
  {
  phi[i * pri[j]] = phi[i] * pri[j];
  break;
  /*
  i % pri[j] == 0换言之,i 之前被 pri[j] 筛过了
  由于 pri 里面质数是从小到大的,所以 i 乘上其他的质数的结果一定也是 pri[j] 的倍数
                 它们都被筛过了,就不需要再筛了,所以这里直接 break掉就好了
  */
  }
  }
  }
 }

注意到筛法求素数的同时也得到了每个数的最小质因子,于是同时求出了欧拉函数。

线性筛法求欧拉函数

注意到在线性筛中,每一个合数都是被最小的质因子筛掉。设 p_1n 的最小质因子,n' = \dfrac{n}{p_1} ,那么线性筛的过程中 n 通过 n' \times p_1? 筛掉。

观察线性筛法的过程,应对 n' mod p_1 分两部分讨论:

  1. n' \bmod p_1 = 0 ????? ,可以说明: n 包含了 n? 的所有质因子。 ??????

    \begin{align*}\varphi(n) &= n \times \prod_{i = 1}^{s} {\dfrac{p_i - 1}{p_i}} \\ &= p_1 \times n' \times \prod_{i=1}^{s}{\dfrac{p_i-1}{p_i}} \\ &= p_1 \times \varphi(n')\end{align*}

  2. n’ \bmod p_i \neq 0??? ,可以说明 i??? 与 p_j??? 互质,根据欧拉函数的性质,有:\begin{align*}\varphi(n) &=\varphi(p_1 \times n') \\&= \varphi(p_1) \times \varphi(n') \\&= \varphi(n') \times (p_1 - 1)\end{align*}???

综上,易得到下面的代码:

 void pre()
 {
     memset(is_prime, 1, sizeof(is_prime));
     int cnt = 0;
     is_prime[1] = 0;
     phi[1] = 1;
     for (int i = 2; i <= MAXN; i++)
    {
         if (is_prime[i])
        {
             prime[++cnt] = i;
             phi[i] = i - 1;
        }
         for (int j = 1; j <= cnt && i * prime[j] <= MAXN; j++)
        {
             is_prime[i * prime[j]] = 0;
             if (i % prime[j])
                 phi[i * prime[j]] = phi[i] * phi[prime[j]];
             else
            {
                 phi[i * prime[j]] = phi[i] * prime[j];
                 break;
            }
        }
    }
 }

例题

  1. HDU2136 - Largest prime factor

  2. HDU4548 - 美素数

  3. HDU6053 - TrickGCD

  4. HDU1286 - 找新朋友

  5. HDU5211 - Mutiple

  6. HDU2710 - Max Factor

  7. HDU6069 - Counting Divisors

  8. CF923A - Primal Sport

  9. CF230B - T-primes

  10. CF45G - Prime Problem

  11. CF449C - Jzzhu and Apples

  12. CF959D - Mahmoud and Ehab and another array construction tasks

  13. POJ3292 -Semi-prime H-numbers

  14. POJ2689 - Prime Distance

  15. POJ2478 - Farey Sequence

  16. POJ2262 - Goldbach's Conjecture

  17. POJ2739 - Sum of Consecutive Prime Numbers

  18. POJ3126 - Prime Path

  19. BZOJ3181 - [Coci2012]BROJ

  20. BZOJ1025 - [SCOI2009]游戏

  21. BZOJ1607 - [Usaco2008 Dec]Patting Heads 轻拍牛头

  22. BZOJ2693 - jzptab