素数的计算: 从试除到筛法(C++实现)

文章目录

什么是素数

素数又称质数,是指大于1的自然数中,除了1和它本身,不能被其它自然数整除的数字。1被定义为非素数。大于1的自然数,如果不是素数则为合数。上古大神欧几里得已经证明了有无穷多个素数存在(欧几里得定理)。

素数定理

引入维基百科的描述

在数论中,素数定理描述素数在自然数中分布的渐进情况,给出随着数字的增大,素数的密度逐渐降低的直觉的形式化描述

素数的个数是有规律可循的,对于正实数x,定义π(x)\pi (x)素数计数函数,即不大于x的素数的个数。
π(x)\pi(x)的估计方法有很多,这里介绍比较简单的一种:

π(x)xlnx\pi(x) \approx \frac{x}{\ln x}

其中lnx\ln x是x的自然对数
维基百科中给出了xx从10到102510^{25}π(x)\pi(x)xlnx\frac{x}{\ln x}数值,可以看出随着xx的增长,π(x)xlnx\frac{\pi(x)}{\frac{x}{\ln x}} 始终小于1.171.17(这里的比值很关键,我们后面需要用来估计用于存储素数的数组的预分配空间的大小)。

常规计算方法(试除法)及实现

基本思想

根据定义,给定一个自然数,先判断它是否是素数。若要计算出所有不大于n的素数,则从2到n,依次判断是否是素数。根据如何判断一个数是否是素数又有以下几种常规方法。

常规方法一

为了判断自然数x是否是素数,我们可以从2开始,分别除以不大于x的每个数,如果存在能整除x的数,则x不是素数。
常规方法一的c++实现

#include <iostream> #include <vector> #include <algorithm> #include <cmath> using namespace std; // 判断n是否为素数 bool is_prime(int n) { if (n < 2) return false; for (int i = 2; i < n; ++i) if (n % i == 0) return false; return true; } // 计算所有不大于n的素数 void get_prime(vector<int>& prime, int n) { for(int i = 2; i <= n; ++i) if(is_prime(i)) // 判断i是否是素数 prime.push_back(i); } int main() { int n = 100000; vector<int> prime; get_prime(prime, n); return 0; }

该算法的时间复杂度为

O(i=2ni2)=O(n2)O(\sum_{i=2}^{n}i^2) =O(n^2)

实际上我们没必要比较这么多次,下面我们介绍方法二。

常规方法二

如果nn为合数,已知n=n×nn=\sqrt{n} \times \sqrt{n},假设n=xyn=xy,如果x>nx>\sqrt{n},那么yy必然小于n\sqrt{n}。即我们只需要判断到n\sqrt{n}即可。
常规方法二的C++实现
实际上只需要根据方法一更改is_prime函数即可,其余代码完全一样。
is_prime函数

bool is_prime(int n) { if (n < 2) return false; /** * 只需要判断到根号n即可,这里必须要使用小于等于符号 * (不能仅仅使用小于符号,例如判断9) **/ for (int i = 2; i * i <= n; ++i) if (n % i == 0) return false; return true; }

该算法的时间复杂度为

O(i=2ni)=O(nn)O(\sum_{i=2}^{n}\sqrt{i}) = O(n\sqrt{n})

实际上我们的算法还有进步的空间,下面介绍常用的两种性能不错的筛法。

埃拉托斯特尼筛法及C++实现

基本思想

埃拉托斯特尼筛法(sieve of Eratosthenes ) 是古希腊数学家埃拉托斯特尼发明的计算素数的方法。对于求解不大于nn的所有素数,我们先找出n\sqrt{n}内的所有素数p1,p2,...,pk1,pkp_1,p_2,...,p_{k-1}, p_k,其中k=nk = \lfloor \sqrt{n}\rfloor,依次剔除pi(1<=i<=k)p_i(1 <= i <=k)的倍数,剩下的所有数都是素数。
维基百科中有一张关于使用埃拉托斯特尼筛法求解素数的GIF图
埃拉托斯特尼筛法求解素数
下面我们来看一个来自维基百科的例子,如果求25内的所有素数
1.列出大于等于2小于等于25的所有数组成的序列:

  • 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25
  1. 去掉素数2的所有倍数:
  • 2 3 5 7 9 11 13 15 17 19 21 23 25
  1. 去掉素数3的所有倍数:
  • 2 3 5 7 11 13 17 19 23 25
  1. 去掉素数5的所有倍数:
  • 2 3 7 11 17 19 23
  1. 因为5=255 = \sqrt{25},所以算法结束。剩下的所有数都是素数。

埃拉托斯特尼筛法的证明

证明之前我们先粗略介绍算术基本定理。算法基本定理的描述是:每个大于1的自然数均可写为素数的积,而且这些素因子按大小排列之后,写法仅有一种方式。然后开始我们的证明(反证法):
假设正整数nn为通过埃拉托斯特尼筛法得到序列中的合数且nn不能被小于等于n\sqrt{n}的素数整除。因为nn是合数,则必然有xy=nxy=n,如果1<x<n1<x<\sqrt{n},那么n<y<n\sqrt{n} < y < n,又因为nn不能被小于等于n\sqrt{n}的素数整除,故xx为合数,根据算术基本定理,合数xx一定写成关于某个素数pp的积,因为1<x<n1 < x < \sqrt{n},则1<p<n1<p<\sqrt{n}。与假设矛盾,证毕。

C++实现

#include <iostream> #include <vector> #include <algorithm> using namespace std; void get_prime(vector<int> &prime, int n) { // 初始化一个布尔数组,用于判断下标i是否是素数 vector<bool> is_prime(n + 1, true); if (n < 2) return; for (int i = 2; i <= n; ++i) { if (is_prime[i]) { prime.push_back(i); // 把素数加入到数组中 // 这里只需要从 i * i 开始筛选即可,证明见下文 for (int j = i * i; j <= n; j += i) is_prime[j] = false; } } } int main() { int n = 100; vector<int> prime; get_prime(prime, n); // 打印计算出的素数 for_each(prime.begin(), prime.end(), [](int &n) { cout << n << " "; }); return 0; }

在上述代码的get_prime函数中

for (int j = i * i; j <= n; j += i)

这一行中的的jj是从i2i^2开始的,也就是说对于当前素数ii是从i2i^2开始筛选的,且剩余序列中从22i21i^2-1的所有数都是素数。证明如下:
如果i=2i=2i2i^2之前的序列为: 2,32,3,都是素数。
如果i>2i>2,因为1<i21<i21<i^2-1 < i^2,所以i1<i21<ii-1<\sqrt{i^2-1} < ii1i-1ii的前一个数,即i21\sqrt{i^2-1}大于i1i-1及之前的所有数。因为ii之前的所有素数都筛选过了,且均小于i21\sqrt{i^2-1},根据埃拉托斯特尼筛法可知,序列中i2i^2之前的所有数都是素数,所以ii只需要从i2i^2开始筛选即可,证毕。
埃拉托斯特尼筛法的时间复杂度为O(nloglogn)O(n\log \log n),证明俺暂时不会(╯-_-)╯╧╧

欧拉筛法及C++实现

基本思想

回顾上文的埃拉托斯特尼筛法,核心思想就是依次剔除素数的倍数。这里我们很容易就能想到此算法存在的缺点:对于两个素数的公倍数,我们可能会进行多次剔除。
于是便出现了欧拉筛法,核心思想便是对于一个没被筛选过的数来说,它只需要被第一个整除它的素数筛掉就行。直接上代码。

C++实现

#include <iostream> #include <vector> #include <algorithm> using namespace std; void get_prime(vector<int> &prime, int n) { if (n < 2) return; // 初始化一个布尔数组,用于记录每个数是否是素数 vector<bool> is_prime(n + 1, true); for (int i = 2; i <= n; ++i) { // 如果是素数则将该数字加入到素数序列中 if (is_prime[i]) prime.push_back(i); // 从第一个素数开始,将它的倍数设置为非素数 for (int j = 0; j < prime.size() && i * prime[j] <= n; ++j) { is_prime[i * prime[j]] = false; /** * 如果 i%prime[j] == 0, 则 i = prime[j] * a * i * prime[j+1] = prime[j] * a * prime[j+1] * prime[j]已经筛过i * prime[j+1]这个数了,不需要用prime[j+1]再筛选一次 **/ if (i % prime[j] == 0) break; } } } int main() { int n = 100; vector<int> prime; get_prime(prime, n); for_each(prime.begin(), prime.end(), [](int &n){ cout << n << " "; }); return 0; }

在上述代码的get_prime函数中

if (i % prime[j] == 0) break;

这里就是欧拉筛法的核心所在啦,对于数字ii来说,如果prime[j]prime[j]能整除ii,则说明i=prime[j]×ai=prime[j] \times a,其中,1<a<i1<a<i,而对于比prime[j]prime[j]大的素数prime[j+1]prime[j+1]i×prime[j+1]=prime[j]×a×prime[j+1]i \times prime[j+1] = prime[j] \times a \times prime[j+1],也就是说i×prime[j+1]i \times prime[j+1]已经被prime[j]prime[j]筛选过了,不需要再筛选了。
欧拉筛法的时间复杂度为O(n)O(n),证明俺暂时也不会(╯-_-)╯╧╧

总结

本文一共介绍了四种计算素数的方法。两种试除法,试除法一时间复杂度为O(n2)O(n^2),试除法二时间复杂度为O(nn)O(n \sqrt{n}),两者的空间复杂度为O(1)O(1);两种筛法,埃拉托斯特尼筛法时间复杂度为O(nloglogn)O(n\log \log n),欧拉筛法时间复杂度为O(n)O(n),两者空间复杂度均为O(n)O(n)
在我们上面的实现中还存在一个问题值得我们思考,那就是如何存储素数,上面的代码中我们是用C++的vector存储素数,vector虽为动态数组,能够动态添加和删除数据,但是在大部分STL的vector的实现版本中,vector都是有预分配空间的,当预分配的空间的容量不够时,便会重新分配更大的空间,并把所有的旧数据复制到新的空间中,所以这里的开销也是难以忽视的(其它语言也存在类似的问题)。所以,我们需要知道我们应该分配多少的空间来存储求得的素数,也就是说我们得估计素数的个数以便预分配存储空间。
而如何估计素数的数量?前面我们介绍的素数定理就排上用场啦,这里就不具体说啦。

参考资料

  1. 质素 - 维基百科,自由的百科全书
  2. 算术基本定理 - 维基百科,自由的百科全书
  3. 素数计数函数 - 维基百科,自由的百科全书
  4. Sieve of Eratosthenes - Wikipedia, the free encyclopedia

本作品采用知识共享署名-非商业性使用-禁止演绎 4.0 国际许可协议进行许可。转载请注明: 作者staneyffer,首发于我的博客,原文链接: https://chengfy.com/post/10


载入评论中....