数数论

素数筛

字数统计: 3.1k字   |   阅读时长: 14分
, , , ,

摘要: 本文介绍素数筛的算法原理和代码模板、以及几道相关题目

【对算法，数学，计算机感兴趣的同学，欢迎关注我哈，阅读更多原创文章】
我的网站：潮汐朝夕的生活实验室
我的公众号：算法题刷刷
我的知乎：潮汐朝夕
我的github：FennelDumplings
我的leetcode：FennelDumplings

各位好，本文详细研究素数筛。首先以力扣第 204 题为模板题学习埃氏筛和线性筛，然后讨论素性测试和分解质因数这两个素数筛的常见应用，分别涉及力扣第 866 题、第 952 题。

素数筛是筛法在素数上的应用，除了素数以外，线性筛还可以用在求积性函数上，例如欧拉函数、莫比乌斯函数。结合算数基本定理，素数筛还可以用于解决约数个数、约数和等问题。

最后例举了一些素数筛的题目，包括 914、1390、LCP14。主要内容如下：

$1 素数筛

首先看以下模板题：

给定整数 n ，返回 所有小于非负整数 n 的质数的数量 。

提示：

1
0 <= n <= 5 * 1e6

示例 1：
输入：n = 10
输出：4
解释：小于 10 的质数一共有 4 个, 它们是 2, 3, 5, 7 。

示例 2：
输入：n = 0
输出：0

示例 3：
输入：n = 1
输出：0

本题就是要求小于 n 有多少个素数。假设有一个函数，输入 n，输出小于 n 的所有素数：

1
vector<int> get_primes(int n);

那么本题的代码就可以写成下面这样：

1
2
3
4
5
6
7
class Solution {
public:
    int countPrimes(int n) {
        vector<int> primes = get_primes(n);
        return primes.size();
    }
};

如何实现 get_primes(n)，这是素数筛要解决的问题。

埃氏筛

如果 x 是合数，那么 x 的倍数也一定是合数。从小到大枚举每个数，然后同时把当前这个数的所有倍数记为合数，结束时未标记的为素数。

代码模板 (C++)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
vector<int> get_primes(int n) {
    if(n < 2) return {};
    vector<bool> vec(n, true);
    vec[0] = false;
    vec[1] = false;
    for(int i = 2; i * i < n; ++i)
    {
        if(vec[i])
        {
            for(int j = i * i; j < n; j += i)
                vec[j] = false;
        }
    }
    vector<int> result;
    for(int i = 0; i < n; ++ i)
    {
        bool flag = vec[i];
        if(flag)
            result.push_back(i);
    }
    return result;
}

时间复杂度证明 $O(\log\log N)$

记 p 为小于 n 的最大素数，内循环的次数如下：

1
2
3
4
5
i == 2: n/2
i == 3: n/3
i == 5: n/5
...
i == p: n/p

总时间复杂度为 $O(\frac{n}{2} + \frac{n}{3} + … + \frac{n}{p})$

下面看 $\frac{1}{2} + \frac{1}{3} + … + \frac{1}{p}$ 为什么等于 $O(\log\log N)$

首先罗列用到的结论：

1. 调和级数 $1 + \frac{1}{2} + \frac{1}{3} + …$

$\lim\limits_{n \to +\infty}\sum\limits_{i=1}\limits^{n}\frac{1}{i}$ 是发散的，但是当 $n \to +\infty$ 时，有近似计算公式 $ln(n) + \gamma$，$\gamma = 0.57721$ 为 欧拉常数。

2. 泰勒级数

上式对 $[-1, 1)$ 的 x 都成立。左右同时取相反数：

3. 欧拉乘积公式

欧拉乘积公式是研究素数分布的基础。欧拉乘积公式的意思是狄利克雷级数可表示为一个与素数相关的无穷乘积。黎曼 $\zeta$ 函数可表示为此无穷乘积的形式。

对任意复数 s， 若 $Re(s) > 1$，则乘积公式如下，其中 n 为正整数，p 为素数。

将 s = 1 代入欧拉乘积公式

两边取对数：

因为 $-1 < \frac{1}{p} < 1$，因此可以对上式泰勒展开: $ln((1 - p^{-1})^{-1}) = \sum\limits_{n=1}\limits^{\infty}\frac{1}{np^{n}}$，而 $p \to \infty$ 时，$\lim\limits_{p \to \infty}\sum\limits_{n=1}\limits^{\infty}\frac{1}{np^{n}}$ 收敛于 $\frac{1}{p}$，因此

将调和级数近似地代入右侧： $\sum\limits_{p}\frac{1}{p} = ln(ln(n) + \gamma)$

最后得到时间复杂度为 $O(\log\log N)$

线性筛(欧拉筛)

埃氏筛的问题在于很多数会重复判断，比如 6 在 i=2 和 i=3 时都判断了一次。

如果能让每个合数都只被标记一次，那么时间复杂度就可以降到 $O(N)$ 了。这就是线性筛。

代码模板 (C++)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
vector<int> get_primes(int n) {
    if(n < 2) return {};
    vector<bool> vec(n, true);
    vec[0] = false;
    vec[1] = false;
    int cnt = 0;
    vector<int> primes;
    for(int i = 2; i < n; ++i)
    {
        if(vec[i])
        {
            ++cnt;
            primes.push_back(i);
        }
        for(int j = 0; j < cnt && i * primes[j] < n; ++j)
        {
            vec[i * primes[j]] = false;
            if(i % primes[j] == 0)
                break;
        }
    }
    return primes;
}

时间复杂度 $O(N)$

关键代码在于 if(i % prime[j] == 0) break;，保证每个合数是被最小质因数筛掉。

因为 primes 数组中的素数是递增的, 当 i 能整除 primes[j]，那么 i * prime[j + 1] 这个合数肯定会被 prime[j] 乘以某个数筛掉。因此，这里直接 break 掉，将 i * prime[j + 1] 及之后的给后面的数去筛。

$2 素数筛应用

素性测试

给定一个正整数，判断该数是否为素数。

1
bool isPrime(int n);

基于试除法

1
2
3
4
5
6
7
8
9
bool isPrime(int n)
{
    if(n == 1)
        return false;
    for(int i = 2; i * i <= n; ++i)
        if(n % i == 0)
            return false;
    return true;
}

基于素数筛

当数据范围较大且需要多次查询的情况下，判断素数可以首先利用筛法筛出所有素数，然后再进行判断：

1
2
3
4
5
6
7
8
bool isPrime(int n, const vector<int>& primes)
{
    int m = primes.size();
    for(int i = 0; i < m && primes[i] * primes[i] <= n; ++i)
        if(n % primes[i] == 0)
            return false;
    return true;
}

866. 回文素数

本题的思路就是构造回文数，然后进行素性测试，题解参考：力扣866-回文素数。

分解质因数

基于试除法

可以扫描 $2\sim\lfloor\sqrt{N}\rfloor$ 的每个 d，若 d 能整除 N，则从 N 中除掉所有的因子 d 并累计除去的 d 的个数。

代码模板 (C++)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
vector<int> p, c;

void divide(int n, vector<int>& p, vector<int>& c)
{
    p.clear();
    c.clear();
    for(int i = 2; i <= sqrt(n); ++i)
    {
        if(n % i == 0)
        {
            // i 是质数
            int cnt = 0;
            while(n % i == 0)
            {
                n /= i;
                ++cnt;
            }
            p.push_back(i);
            c.push_back(cnt);
        }
    }
    if(n > 1)
    {
        // n 是质数
        p.push_back(n);
        c.push_back(1);
    }
}

一个合数，其因子一定在扫描到该合数之前就被筛掉了，即在试除的过程中，能整除 N 的一定是质数。时间复杂度 $O(\sqrt{N})$。

在题目 【分治】A的B次方的所有约数之和 中，我们实践了上述分解质因数的代码模板。

基于素数筛

如果已经通过素数筛得到了 $2\sim\lfloor\sqrt{N}\rfloor$ 的素数集合 primes。可以用以下代码，返回的 result 是质因数。

代码模板 (C++)

  • 返回的 result 无重复
1
2
3
4
5
6
7
8
9
10
11
12
13
vector<int> get_prime_factor(int num, const vector<int>& primes) {
    // 大于 1 的质因数
   if(num <= 0) return {};
   if(num == 1) return {};
   vector<int> result;
   for(int prime: primes)
   {
       if(prime > num) break;
       if(num % prime == 0)
           result.push_back(prime);
   }
   return result;
}
  • 返回的 result 有重复
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
vector<int> get_prime_factor(int num, const vector<int>& primes) {
    // 大于 1 的质因数
   if(num <= 0) return {};
   if(num == 1) return {};
   vector<int> result;
   int m = primes.size();
   int i = 0;
   while(i < m && primes[i] <= num)
   {
       if(num % primes[i] == 0)
       {
           result.push_back(primes[i]);
           num /= primes[i];
       }
       else
           ++i;
   }
   return result;
}

Pollard-rho 算法是比试除法更高效的算法，属于随机算法。

952. 按公因数计算最大组件大小

题解参考：力扣952-按公因数计算最大组件大小。

$3 筛法应用

线性筛基本可以求出所有积性函数，方法会有相同，外层循环都是从2开始的。

积性函数: 对于所有互质的整数 a 和 b 有性质 $f(ab)=f(a)f(b)$ 的数论函数。

筛法求欧拉函数

欧拉函数是 <= n 的正整数中与n互质的数的数目, 例如 $\phi(1)=1$，$\phi(8) = 4$ (1,3,5,7与8互质)。

其中 $p_{1}, p_{2}, …, p_{k}$ 为 n 的 k 个质因数。推导需要用到容斥原理，见 组合数学4。

性质：

  • (1) $\phi(n)$ 是积性函数：$当 n 和 m 互质时，\phi(mn) = \phi(n)\phi(n)$
  • (2) n 为质数则 $\phi(n) = n - 1$
  • (3) n 为奇数则 $\phi(2n) = \phi(n)$
  • (4) 若 a 为质数, $b \equiv a \mod 0$, 则 $\phi(ab) = a\phi(b)$
  • (5) n 为 $p^{k}$ 则 $\phi(n) = p^{k} - p^{k-1}$
  • (6) $\sum\limits_{d|n}\phi(d) = n$

欧拉函数是积性函数，因此可以用筛法计算出某个范围内所有数的欧拉函数值，在素数线性筛的基础上稍微修改。

代码模板 (C++)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
vector<int> get_phi(int n)
{
    if(n < 2) return {};
    vector<bool> vec(n, true);
    vec[0] = false;
    vec[1] = false;
    int cnt = 0;
    vector<int> primes;
    vector<int> phi(n, -1);
    phi[1] = 1;
    for(int i = 2; i < n; ++i)
    {
        if(vec[i])
        {
            ++cnt;
            phi[i] = i - 1; // 性质2
            primes.push_back(i);
        }
        for(int j = 0; j < cnt && i * primes[j] < n; ++j)
        {
            vec[i * primes[j]] = false;
            if(i % primes[j] == 0)
            {
                phi[i * primes[j]] = phi[i] * primes[j]; // 性质4
                break;
            }
            else
                phi[i * primes[j]] = phi[i] * (primes[j] - 1); // 性质1: 积性函数
        }
    }
    return phi;
}

筛法求莫比乌斯函数

莫比乌斯函数:

从以上定义可以知道，莫比乌斯函数与非平方数的质因子数目有关。以上定义的通俗表达：

  • $\mu(1) = 1$
  • 当 n 存在平方因子时，$\mu(n) = 0$
  • 当 n 是素数或奇数个不同素数之积时，$\mu(n) = -1$
  • 当 n 是偶数个不同素数之积时，$\mu(n) = 1$

$\mu(n)$ 是积性函数。可以用筛法求某个范围内所有的莫比乌斯函数值。

代码模板 (C++)

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
vector<int> get_mu(int n)
{
    if(n < 2) return {};
    vector<bool> vec(n, true);
    vec[0] = false;
    vec[1] = false;
    int cnt = 0;
    vector<int> primes;
    vector<int> mu(n, -2);
    mu[1] = 1;
    for(int i = 2; i < n; ++i)
    {
        if(vec[i]) // i 是质数
        {
            ++cnt;
            primes.push_back(i);
            mu[i] = -1;
        }
        for(j = 0; j < cnt && i * primes[j] < n; ++j)
        {
            vec[i * primes[j]] = false;
            if(i % primes[j] == 0)
            {
                mu[i * primes[j]] = 0;
                break;
            }
            else
                mu[i * primes[j]] = -mu[i];
        }
    }
}

筛法求约数个数

$n = p_{1}^{a_{1}}p_{2}^{a_{2}}…p_{k}^{a_{k}}$, 约数个数 $fac(n) = (1 + a_{1})(1 + a_{2})…(1 + a_{k})$

筛法求约数和

$n = p_{1}^{a_{1}}p_{2}^{a_{2}}…p_{k}^{a_{k}}$, 约数和 $sumfac(n) = \sum\limits_{i=0}\limits^{a_{1}}p_{1}^{i}\sum\limits_{i=0}\limits^{a_{2}}p_{2}^{i}…\sum\limits_{i=0}\limits^{a_{k}}p_{k}^{i}$

$4 题目

914. 卡牌分组

题解参考文章：最大公约数算法以及C++17组件

1390. 四因数

题解参考文章：力扣1390-四因数

LCP 14. 切分数组

题解参考文章：力扣LCP14-切分数组

Share