基于 MPI 实现埃拉托斯特尼筛法及性能优化

1 项目背景
1.1 素数
素数，又叫质数，是指一个大于 1 的自然数，且除了 1 和它本身外，不能被其他自
然数整除的数。换句话说，就是该数除了 1 和它本身以外，不再有其他的因数。
在初等数学中有一个基本定理，任意一个大于 1 的自然数，要么本身就是素数，要
么可以分解为几个素数之积，这种分解本身，具有唯一性。基于此，现在多将素数数用
于密码学上，而其解密的过程，实际上就是一个寻找素数的过程。
素数的常见求解方法有：
• 筛法方法 (Sieve method)：包括埃拉托斯特尼筛法和其他素数筛选法。这些方法
通过筛选和排除合数来确定素数。
• 费马测试 (Fermat’s test): 基于费马小定理，随机选择数 a 进行测试，如果不符合
定理，该数不是素数。
• 米勒-拉宾素性测试 (Miller-Rabin primality test)：一种概率性的素数测试方法，通
过选择随机数 a 进行多次测试，如果不符合定理，该数不是素数。
• 素数判定法 (Primality testing methods)：包括传统的试除法和更高效的方法，如
费马素性检验、米勒-拉宾素性检验和 Solovay-Strassen 素性检验等。
• 梅森素数检查（Mersenne primality test）：用于确定梅森数（形如 2 p - 1）是否为
素数的方法。这种方法结合了特殊的数学性质和梅森素数的特定形式来进行检查。
这些方法在实际中被广泛应用，每种方法都有其优点和适用范围。选择合适的方法
取决于特定的求解问题和所需的准确性要求。
素数的应用主要在以下 4 个方面：
• 被利用在密码学上，所谓的公钥，就是将想要传递的信息在编码时加入质数，编
码之后传送给收信人，任何人收到信息后，若没有此收信人所拥有的密钥，则解
密的过程中（实为寻找质数的过程），将会因为找分解质因数过久，而失去时效性。
• 在汽车变速箱齿轮的设计上，把相邻的两个大小齿轮的齿数设计成质数，以增加
两个相同的齿相遇啮合次数的最小公倍数，可增强耐用度减少故障。
• 以质数形式无规律变化的导弹和鱼雷可以使敌人不易拦截。
• 多数生物的生命周期也是质数（单位为年），这样可以最大程度地减少碰见天敌的
机会。
1.2 云计算平台
云计算平台也称为云平台，是指基于硬件资源和软件资源的服务，提供计算、网络
和存储能力。云计算平台可以划分为 3 类：以数据存储为主的存储型云平台，以数据处
理为主的计算型云平台以及计算和数据存储处理兼顾的综合云计算平台。云计算平台
具有降低计算机成本、优化性能、降低软件成本、提高数据可靠性等诸多优点。利用华
为云提供的弹性云服务器平台配置 MPI 环境，可以在云端实现项目的构建和调试、运
行，减少对于本地硬件的需求，提供了更高的灵活性。
1.3 MPI
MPI(message passing interface) 消息传递接口是一个为并行计算设计的通信协议，
支持点对点和广播。MPI 是一个信息传递应用程序接口，包括协议和语义说明，他们指
明其如何在各种实现中发挥其特性。MPI 的目标是高性能，大规模性和可移植性。MPI
在今天仍为高性能计算的主要模型，是一种基于信息传递的并行编程技术；消息传递
接口是一种编程接口标准，而不是一种具体的编程语言。简而言之，MPI 标准定义了
一组具有可移植性的编程接口，在当前的并行机上，比较流行的并行编程环境，可以分
为三类：消息传递，共享存储，和数据并行。与这三类对应的典型代表分别为：（MPI，
PVM），OpenMP，和 HPF（High Performance Fortran）。相比其他的并行编程环境，
MPI 具有很多优点，具有很好的可移植性，很好的扩展性与可伸缩性，比其他消息传递
系统好用，有完备的异步通信功能，有精确的定义。
2 原理分析
埃拉托斯特尼筛法（sieve of Eratosthenes），简称埃氏筛，也称素数筛，是简单且
历史悠久的筛法，用来找出一定范围内所有素数。
在寻找整数Ｎ以内的素数时，古希腊数学家埃拉托斯特尼采用了一种与众不同的
方法：先将２－Ｎ的各数写在纸上：
在２的上面画一个圆圈，然后划去２的其他倍数；第一个既未画圈又没有被划去的
数是３，将它画圈，再划去３的其他倍数；现在既未画圈又没有被划去的第一个数是５，
将它画圈，并划去５的其他倍数……依此类推，一直到所有小于或等于Ｎ的各数都画了
圈或划去为止。这时，画了圈的以及未划去的那些数正好就是小于Ｎ的素数。
以 25 为例，详细列出算法如下：
• 列出 2 以后所有数：
2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25
• 标记第一个质数 2
• 标记 2 的倍数
• 如果最大数不大于最后一个标出的素数的平方，那么剩下的所有的数都是质数，
否则回到第二步
本例中，25 大于 2 的平方，返回第二步；2 之后第一个质数是 3，标记 3 的倍数；
得到的质数是 2、3，25 仍大于 3 的平方，再次返回第二步；3 之后第一个质数是 5，标
记 5 的倍数；得到的质数是 2、3、5，25 是 5 的平方，筛选完毕。
去掉标记的数，25 以内的质数是
2，3，5，7，11，13，17，19，23

给出其 C 语言串行实现：
`
int prime[100005];

bool is_prime[1000005];

int eratosthenes(int n) {
int p = 0;
for (int i = 0; i <= n; i++) {
is_prime[i] = true;
}
is_prime[0] = is_prime[1] = 0;
for (int i = 2; i <= n; i++) {
if (is_prime[i]) {
prime[p++] = i;
if (1ll * i * i <= n) {
for (int j = i * i; j <= n; j += i) {
is_prime[j] = 0;
}
}
}
}
return p;
}

`
下面使用并行计算加快这一筛法将数组分为 p 个连续的块, 每个块大小基本相等。
平衡负载要给每个进程分配 ⌊ n
p ⌋ 或 ⌈ n
p ⌉ 个元素, 我们考虑下面这种不同的实现方法：
进程 i 控制的第一个元素是 ⌊in/p⌋，最后一个元素是 ⌊(i + 1) n/p⌋ − 1，对于特定
元素 j，控制他的进程是 ⌊(p (j + 1) − 1) /n⌋。
我们将数组 {2, 3, · · · , n} 分配给 p 个进程，进程 id 分到的数据块为
{2 + ⌊id ∗ (n − 1) /p⌋, 3 + ⌊id ∗ (n − 1) /p⌋, · · · , 1 + ⌊(id + 1) ∗ (n − 1) /p⌋}
令
low_value = 2 + ⌊id ∗ (n − 1) /p⌋
high_value = 1 + ⌊(id + 1) ∗ (n − 1) /p⌋
size = ⌊(id + 1) ∗ (n − 1) /p⌋ − ⌊id ∗ (n − 1) /p⌋
52 原理分析
0 进程分配到的数据块大小 proc0_size = ⌊(n − 1) /p⌋
我们用进程 0 来储存步骤中用于筛选的 k（即 2 到 √ n 的质数），所以程序运行的
前提必须要求
2 + proc0_size >= (int)sqrt((doubt)n)
对每一个进程都要提供一个 marked[size] 这样的数组，prime 保存当前用于筛选的质数，
first 表示进程 id 中第一个要被筛掉的数对应的 marked 数组的下标。marked[size] 数组
代表 {2 + ⌊id ∗ (n − 1) /p⌋, 3 + ⌊id ∗ (n − 1) /p⌋, · · · , 1 + ⌊(id + 1) ∗ (n − 1) /p⌋} 被标记
的情况。
index 用于步骤中找到比 prime 大的未被标记的数中最小的那个数，index 为 0 进
程专属。核心部分的程序为：

`

if(!id) index=0;//0进程专属

prime=2;//最开始用2去筛选

do{

    //找到第一个被prime筛掉的数
    
    if(prime*prime>low_value) first=prime*prime-low_value;
    
    else if(!(low_value%prime)) first=0;
    
    else first=prime-(low_value%prime);
    
    //筛数，筛掉的数对应marked数组相应位置赋值为1，做标记
    
    for(i=first;i<size;i+=prime) marked[i]=1;
    
    //在0进程中找到比prime大的未被标记的数中最小的那个数，用这个质数做新的prime去筛选。
    
    if(!id){
    
        while(marked[++index]);
        
        prime=index+2;//marked[0]对应自然数2
        
    }
    
    //把0进程找到的新的prime更新其他进程的prime的值
    
    MPI_Bcast(&prime,1,MPI_INT,0,MPI_COMM_WORLD);
    
}while(prime*prime<=n);

结合 MPI 并行运算分析可能的几种并行化改造方案：
• 特定线程剔除特定数：优点：直观；缺点：2 的倍数最多，3 的倍数次多，前几个
进程几乎决定了总的运算时间，成为瓶颈。
• 按照可用进程数分段：将待筛选的数字分为不同的段，每个段包含连续的一部分
数字。首先，每个进程独立地筛选其分配到的段中的素数。接下来，进程间进行通
信，共享各个进程的筛选结果，以便确定进一步筛选的范围。重复上述步骤，直
到完成筛选过程。
根据流程进行时间复杂度分析：
63 技术路线
χ : 执行二元操作所需的时间。
λ : 经由通道将一个整数传到另一个个通道所需的时间。
素数定理：π (x) ∼ ln x x，其中 π (x) 表示不超过 x 的素数个数。
2 到 √ n 内，找到一个质数需要一次广播，所以通信的预期时间为 (
√ n/ ln √ n) λ log
2
p，
对各质数 p，需要筛 n/p 次，一共需要 ∑ p n/p，又因为素数倒数和
∑ p 1/p ∼ O (ln ln n
)
所以 p 个进程的筛选的预期时间粗略估计为
χ
(n ln ln n) /p，总预期时间为
χ (n ln ln n) /p + (√ n/ ln √ n ) λ log2 p
由此可以想到代码优化方案，如删除偶数，这样筛选的预期时间约为 χ (n ln ln n) /2p；
消除广播，使用串行算法求出 2 到 √ n 的质数，最后通信归总结果，只需要 λ⌈log2 p⌉
的通信时间，一共需要 λ⌈log2 p⌉ + χ [(n ln ln n) /2p + √ n ln ln √ n]。
3 技术路线
3.1 环境搭建
3.1.1 弹性云服务器
弹性云服务器（Elastic Cloud Server，ECS）是由 CPU、内存、操作系统、云硬盘
组成的基础的计算组件。弹性云服务器共分为通用计算型、内存密集型、高性能计算型、
计算加速型等十数种规格实例类型，满足各种场景的上云需求。弹性云服务器计费模式
灵活，数据可靠，提供多维度安全防护，弹性易用，运维高效，实时云端监控，负载均
衡。我们的整个项目程序均通过华为云弹性云服务器运行。
3.1.2 环境配置步骤

1. 购买云服务器
2. 配置节点互信
3. 安装 MPI
   关键节点截图如下图所示。
   3.2 程序编译与脚本运行
   3.2.1 程序编译指令
   1 #!/bin/bash
   2 compile_command() {
   3
   ver=$1
   4
   cpp_file=2 83.3 版本优化介绍 3 技术路线 5 cho "Compiling ver_verfrom cpp_file" 6 mpicxx -g -Wall -o ver_ver cpp_file 7 echo "---------------------" 8 } 9 # 循环执行编译命令 10 for ver in {1..3}; do 11 cpp_file="ver_ver.cpp"
   12
   compile_command ver
   cpp_file
   13 done
   3.2.2 脚本批量运行
   1 #!/bin/bash
   2 execute_command() {
   3
   ver=$1
   4
   np=$2
   5
   n=$3
   6
   7
   echo “Executing ver_verwith−np
   np and n = $n”
   8
   mpiexec -f config -np np ./ver_
   ver $n
   9
   echo “---------------------”
   10 }
   11 # 循环执行不同的参数组合
   12 for ver in {1…3}; do
   13
   for np in {1…8}; do
   14
   for n in 10000 1000000 100000000 1000000000; do
   15
   execute_command ver
   np $n
   16
   done
   17
   done
   18 done
   3.3 版本优化介绍
   • Ver1.0：定义一个数组 marked, 每一个元素的下标对应一个整数，它的值表示这
   个整数是否为素数, 值为 1 是素数，值为 0 不是素数。先假定所有的数都是素数，
   将 marked 数组置 0。选定第一个整数２，开始依次标记 2 的倍数，一直标记到
   最后一个数为止。接下来选定接下来第一个未标记的数，它一定是素数，在使用
   广播的形式通知各进程筛选出这个素数的倍数。这样循环到最后，所有进程中未
   93.4 运行结果分析
   3 技术路线
   标记的数之和就是 1-n 中的所有素数了。
   • Ver2.0: 利用已知除 2 以外的所有偶数都不是素数的常识，可以将待筛选数字总量
   减半，从而提高筛选效率。关键代码在于数组减半，找到新的索引映射，以及首
   个倍数（非素数）的位置。
   • Ver3.0: 初始的代码是通过进程 0 广播下一个筛选倍数的素数。进程之间需要通过
   MPI_Bcast 函数进行通信。通信会有一定开销，特别是在分布式计算机架构上，
   因此我们让每个进程都各自找出它们的前 sqrt(n) 个数中的素数，在通过这些素
   数筛选剩下的素数，这样一来进程之间就不需要每个循环广播素数了，性能得到
   提高。
   • Ver4.0: 当计算规模较大时使用合理设置缓存大小优化思路原理本优化方法需每
   个进程都各计算出前 sqrt(n) 个数中的素数。外存储读取耗时远大于直接操作内
   存，设置缓存优化的思路在于每次处理缓存大小的数组，之前我们已经将 n 内分
   成大小约为 n/p 的块给每个进程处理，然后在在每个进程中将 n/p 大小块按照缓
   存大小进行分块，另外根据测试的时候 n 大小进行选择 cache 的级别，比如测试
   的是亿级的数据，远超过 cache 的大小，所以直接对 L3 级别的 cache 进行分块，
   当然选择 L3 并不一定是最优策略，需要多次实证才能知道。这里需要使用 linux
   查看所使用云服务器 CPU 的缓存大小进行合理设置，命令如下:
   1 lscpu | grep “cache”
   代码均已放入附录部分。
   3.4 运行结果分析
   3.4.1 结果展示
   如图五所示, 列出了在 1E5、1E7、1E9、1E10 四种数量级下三种不同算法的耗时
   对比，并在 1E10 中加入了设置缓存后的耗时对比。
   3.4.2 分析
   运行耗时随着启用的线程数增加而减少是典型的并行计算的预期行为。然而，当计
   算难度较小时，发现运行耗时与线程数没有明显的改善，可能是由以下原因引起的：
   • 任务划分不均匀：当计算难度较小时，任务的计算量可能相对较小。如果任务在
   划分时不均匀，导致某些线程的计算负载较轻，而其他线程的计算负载较重，那
   么耗时的波动性可能会增大，并且无法充分利用所有线程的计算能力。
   同步开销：在并行计算中，线程之间需要进行通信和同步。当计算任务较小时，线
   程之间的通信和同步开销可能会显得相对较高，导致运行耗时没有明显的改善。
   这是因为通信和同步操作本身需要时间，并且随着线程数的增加，这些开销可能
   会更加显著。
   • 硬件限制：在某些情况下，硬件资源可能成为瓶颈，限制了并行计算的性能提升。
   例如，当计算任务较小时，线程间的通信和同步可能会受到网络带宽或内存带宽
   的限制，从而无法实现理想的加速比。当计算难度较大时，耗时随着线程数增大
   而明显减少的情况可能是因为任务更加密集，线程的计算负载更加均衡，同时并
   行计算能够更好地利用硬件资源和并行计算的优势。此时，通信和同步的开销可
   能相对较小，而计算部分的耗时占据了主要的时间。
   为了更好地利用并行计算的优势，可以尝试以下方法：
   • 确保任务划分均匀，使得每个线程的计算负载相近
   • 使用更高效的通信和同步机制，减少开销。例如，可以尝试减少通信次数、使用
   非阻塞通信等技术。
   • 优化算法和计算过程，减少计算量，提高计算效率。
   114 总结展望
   • 考虑硬件资源的限制，合理设置线程数，避免超过硬件能力的范围。
   通过优化任务划分、通信和同步机制以及算法设计，可以提高并行计算的性能，并
   实现更好的加速比。
   4 总结展望
   本次课程项目通过对埃拉托斯特尼筛法的原理分析，将其与 MPI 并行运算结合，
   大大降低了运算时间和时间复杂度。
   在此基础上，我们又将时间复杂度具体化，分析了各部分对于时间复杂度的影响，
   由此提出了删除偶数和消除广播两种降低时间复杂度的方法，在 1E5、1E7、1E9、1E10
   四个量级的数据上基本都能降低一半的运行时间。
   根据最终运行时间图像可以发现，在一定规模的数据量下，随着进程数量的增加，
   运行时间是逐渐减小的，同时在 1E9 和 1E10 规模的图像中发现，去除偶数的算法耗费
   的时间略高于优化通信的时间，说明任务在线程中通信花费的时间较多，在 1E5 规模
   下时间随进程增大也能说明这一结论。