我们的求质数终于进行到万亿了. 不过在这之前, 我们求质数一直都是用的单核单线程程序, 那么现在, 我们用一下多线程呗.

本篇文章主要介绍用C语言利用多线程计算一万亿( 10^{12} )以前的所有质数的算法和代码, 储存质数的算法和代码见下一篇文章.

算法是我们以前求 100 亿质数表中用到的埃拉托色尼筛[1] ,

这种筛法的原理是先将 2\sim n 的所有数字列出来.

\begin{matrix}2&3&4&5&6&7&8&9&10&11&12&13&14&15&16&17\end{matrix}

从小往大阅读. 第一个数字是 2 , 没有被筛掉. 所以第一个质数是 2 . 然后我们把所有 2 的倍数筛掉.

\begin{matrix}\color{blue}{2}&3&\cancel4&5&\cancel6&7&\cancel8&9&\cancel{10}&11&\cancel{12}&13&\cancel{14}&15&\cancel{16}&17\end{matrix}

第二个数字是 3 , 没有被筛掉, 所以第二个质数是 3 , 然后我们把所有 3 的倍数筛掉.

\begin{matrix}\color{blue}2&\color{blue}3&\cancel4&5&\cancel6&7&\cancel8&\cancel{9}&\cancel{10}&11&\cancel{12}&13&\cancel{14}&\cancel{15}&\cancel{16}&17\end{matrix}

第三数字是 4 , 被筛掉了, 我们读取下一个, 5 , 没有被筛掉. 所以第三个质数是 5 .

但这时候, 5>\sqrt{17} , 所以不需要再筛了. 剩下的都是质数. 因此我们输出所有没有被筛掉的数.

\begin{matrix}2&3&5&7&11&13&17\end{matrix}

这些就是 [2,17] 的所有质数.

所以, 如果要用单核单线程的程序求 10^{12} 以前的所有质数, 就只需要先求出 10^{6} 的所有质数(下文简称 质因子),

再让每个质因子划去他所有的在 10^{12} 以内的倍数.

单线程情况如此. 多线程我们要怎么做呢? 很简单, 埃拉托色尼筛是非常适合并行计算的. 我们让每一个线程处理一个质因子就好了. 当然, 一共有 78498 个质因子, 我们一般的计算机可没有这么多核心数, 我们只能创建一点点线程, 当一个线程处理结束后, 去拿下一个质因子, 当然, 其他线程拿的, 他就不需要再拿了.

总结一下:

前面提到, 埃拉托色尼筛的核心算法是筛掉合数, 那么我们必须用一定的内存空间来保存一个数是否被筛去.

最简单的, 就是筛去一个数的时候, 就赋值 0 给他.

但是, 一万亿, 10^{12} 已经大于 unsigned 能最大表示的数 2^{32}=4294967296 了.

所以我们只能用 long long unsigned int 来储存质数.

但是, 一个 long long unsigned int 长达 8 字节. \begin{align} \frac{8\times10^{12}}{2^{30}}=7450.58 \end{align} . 也就是说, 我们需要用 7450.58\,\text{GiB} 的内存来保存这些数的信息.

这肯定在现阶段是不太现实的.

因此, 我们用一下位运算, 只用一位来表示一个数是否被划去. 这样一来, 我们就只需要

\begin{align} \frac{10^{12}}{8\times2^{30}}=116.4153\,\text{GiB} \end{align} 的内存了.

但是, 好巧不巧, 就是没有这么多内存. 所以我们再去掉所有的偶数, 反正偶数除了 2 以外也肯定不是质数. 这样一来, 我们就只需要 \begin{align} \frac{10^{12}}{2\times8\times2^{30}}=58.2077\,\text{GiB} \end{align} 的内存了.

为了方便管理和用多线程, 我们将这个 \begin{align} 59\,\text{GiB} \end{align} 的空间用二维数组保存.

(注意, 该程序仅在Linux或类Unix操作系统上编译和运行)

我们再创建一个长度为 59 的 long long unsigned int 指针数组, 这个数组叫做 primArray , 每个元素指向一个长度为 1\,\text{GiB} 的 long long unsigned int 数组, 并且用 indexArray 表明这是第几个数组.

为每一个 indexArray 数组申请内存.

注意, 因为在这个程序中, 线程函数要传的参数实在是太多了, 所以那些多线程函数要用的数据, 我都定义为了全局变量.

在main函数中创建一个 pthread_t 数组, 用来储存即将创建的线程.

使用 pthread_create 函数创建线程. 不传参(都被定义为全局变量在函数中直接调用了), 所以最后一个参数为NULL.

需要每个线程拿到自己的编号之后, 才能继续. (其实这里最好用互斥量, 但这个代码是我用后面的代码搬过来的, 所以懒得改了)

函数 threadMalloc 在59个线程上同时执行, 向操作系统申请所需的内存空间.

这里申请内存空间的代码, 如果是单线程执行, 需要大概一分钟, 59个线程同时执行的话, 只需要大概2秒.

等待 threadMalloc 函数都执行完毕才进行以后得操作.

下一步, 是不是让线程获取自己将要处理的质因子呀? 是的, 但是在这之前, 我们先用单线程求出 10^{6} 之前的所有质数, 反正这个速度很快, 没必要用多线程.

这样做的目的是, 防止有些本来是合数的数, 因为没有来得及被其他线程划去, 而被当成了质数被另一个线程取出来. 这种情况不会导致结果错误, 但会消耗无意义的时间.

我们用一个简单的函数来处理这个需求.

下一步, 就是每个线程拿取自己的质因子, 并且划去质因子的倍数的函数了.

(如果不知道这个算法具体是怎么实现的, 请看参考链接1)

其中, 最重要的是 getOut 函数, 用来特定地划去一个数.

求出这个数所在的 indexArray 也就是处于 primArray 的哪一个 1\,\text{GiB} 的大小的数组中.

求出这个数处于 primArray[indexArray] 所指向的1\,\text{GiB} 的大小的数组中的哪一个 long long unsigned int 中.

求出在 long long unsigned int 的哪一位中.

(一个 long long unsigned int 有 8 个字节, 所以可以保存 64 个二进制位信息).

之后, 再把这一位标为 0 , 就代表他已经被划去了.

对了, 对于每个线程来说, 都要拿取自己的质因子, 所以我们有一个函数叫做 getNextPrime

取到当前的质因子之后, 要保存到全局变量 readFromBitsStruct->thisPrime 中, 这样其他线程在取质因子之前, 就知道上个线程取走了哪一个, 他应该取下一个.

写好这个函数之后, 就在main函数中创建这些线程, 再等待他们都运行结束.

这样就好了? 我们来思考一下有什么问题?

当两个线程同时取质因子, 有可能取到的是同一个质因子. 这个不会对结果造成错误, 但是会浪费时间.

所以我们在取质因子的时候, 加一个互斥量 mutex1, 取完之后, 让这个互斥量+1, 允许其他线程取质因子. 如果互斥量等于1, 就代表线程可以进入这个临界区, 如果等于0 , 就得等着其他线程退出临界区他再进去.

当然, 在用这个互斥量之前, 得初始化这个互斥量.

同时注意, 互斥量只能是全局变量

这样就结束了吗?

其实现在已经可以运行了.

但运行结束之后会发现质数的个数多了大概 \begin{align} \frac{1}{1000} \end{align} , 为什么呢?

我们思考一下, 在 getOut 函数中, 标记一个质数不是质数的语句, 究竟在计算机里面是怎么执行的.

每划去一个数之前, 都要算出管理这个数的是哪一个信号量.

之后, 在处理之前, 都得先查询一下信号量是不是为 1 . 如果不是 1 就等一下. 每一个线程处理之后, 都要 sem_post 一下信号量, 让信号量的值变回 1 , 允许下一个进程进入.

另外, 为了防止在刚刚开始的时候, 大量线程都在等待同一个信号量. 我们在创建线程之后阻塞主线程 1000 微秒, 之后才能创建下一个线程. 让刚刚创建的线程先跑一段时间, 跑到前面去.

这样, 结果就正确了.

在运算结束后, 还需要读出质数, 所以我们写一个叫做 threadReadFromBits 的函数

当然, 为了方便将这些质数写成文件, 我们得先释放一部分内存. 不过在释放之前, 最好销毁掉.

之后, 我们就按位判断每个数到底是不是质数, 如果是质数, 还要判断是不是每个 indexArray 的第一个质数.

如果是质数, priemCount_inner 要加加, 最后, 要写到 RFBS_Thread 这个数据结构里面去, 这个数据是存储质数个数和最大的质数的.

经过运算, 在一万亿( 10^{12} )之前, 一共, 37607912018 个质数,

第一个质数是 2 , 最后一个质数是 999999999989 , 和已经验证的数值是一样的.

一共需要大概一个半小时的时间. 时间主要花在线程等待信号量上面, 因为如果不加信号量, 虽然质数的个数会算出来多千分之一, 但是只需要大概半个小时.

(双路Xeon E5 2667V4, 我内存速度有点慢, 如果内存速度快的话会更快)

我们算一下多线程倍率?

32 线程要跑 4626 秒, 多线程要跑 63135 秒, 所以多线程倍率

\begin{align} \frac{63135}{4626}=13.6479 \end{align}

本程序已开源, 开源地址:

第二篇文章: 介绍储存质数的算法和代码

运行本程序的机器:

写在后面:

"质数" 系列文章, 未来还会继续推进吗?

应该在未来很长一段时间, 都不会再推进质数计算了, 一来一万亿以前的所有质数, 在很多场景下, 是完全够用的,

二来, 我们从试除法, 一直到筛法, 一直到多线程, 其实现在计算质数的方法, 没有多少比这些更高效的了.

最重要的是, 如果下一次要计算 10 万亿之前的所有质数, 按照现在的方法, 至少需要 930\,\text{GiB} 内存.

减少内存占用的方法有, 我们可以处理一段将内存写回磁盘, 但那样效率就比较低了, 而且储存空间的占用也是需要考虑的.

最后, 非常感谢我的读者, 如果不是你们一直以来支持我, 我应该不会把质数系列推进到现在.

也非常感谢一路走来, 给我提供建议的读者, 我从你们的建议中学到了非常多的方法和技巧.

非常感谢你们, 我写的文章能被你们看到, 或许, 这就是 yuanfen 吧!

封面图片

クリスマスおうちデート