MPI并行埃拉托斯特尼筛法：从串行到分布式的性能优化实战

简介：这份资源面向学习并行计算、MPI编程及C++课程设计的高校学生与开发者，围绕埃拉托斯特尼筛法这一经典素数求解算法，提供从串行实现到MPI多进程并行优化的完整代码方案，帮助读者理解素数筛选原理与并行任务划分、通信开销等核心问题。压缩包共356个文件，约1.15MB，以276个txt说明与测试数据、28个cpp源码、9个cmake构建脚本为主，另含bin与exe可执行文件、csv结果数据、make与makefile构建配置、docx文档及少量C、Python、Shell辅助脚本，覆盖编译、运行、验证全流程。资源已有380人学习下载，适合作为课程设计参考或并行编程入门练手项目。读者可获取可直接编译运行的MPI筛法实现、构建配置与结果数据，对照串行与并行版本分析加速比与性能瓶颈，并借助文档与脚本快速复现实验、排查编译与运行问题。

1. 从串行到并行：MPI 埃拉托斯特尼筛法到底能压榨出多少性能

单机筛 10 亿以内的素数，串行 C 程序大概要跑十几秒，内存占用 1 GB 左右。这个数字放在教学场景里没问题，但一旦数据规模上到百亿级，或者需要在集群上做分布式素数统计，串行版本就彻底不够看了。MPI 埃拉托斯特尼筛法要解决的核心问题就是：把筛法的标记过程拆到多个进程上，让每个进程只负责一段区间，通过通信协调边界素数的广播，最终把整体耗时压到可接受的范围。这不是简单的“把循环拆开”就能搞定的事——筛法的本质是全局依赖，每个合数都要被它最小的素因子标记，而最小素因子可能来自任意进程负责的区间。所以 MPI 版本的关键在于：如何划分区间、如何收集和广播素数、如何避免重复标记和通信爆炸。适合有 C/C++ 基础、了解 MPI 基本点对点通信和集合通信、想从串行算法过渡到分布式算法的工程师。如果你正在做大规模数论计算、密码学预处理或者只是想搞懂 MPI 在规则计算场景下的性能边界，这个方向值得花时间。

2. 筛法并行化的理论边界：为什么不能简单对半砍

2.1 埃拉托斯特尼筛法的串行依赖到底卡在哪

串行筛法的逻辑很直白：从 2 开始，把每个素数的倍数全部标记为合数，剩下的就是素数。伪代码大概长这样：

// 串行筛法核心逻辑 for (int i = 2; i * i <= N; i++) { if (!composite[i]) { for (int j = i * i; j <= N; j += i) { composite[j] = 1; } } }

这段代码的依赖关系是：标记j的时候，必须已经知道i是素数。而i是不是素数，又依赖于之前所有比i小的素数有没有标记过它。这是一个典型的前缀依赖链。如果你把[2, N]切成两半，进程 0 负责[2, N/2]，进程 1 负责[N/2+1, N]，进程 1 在标记的时候需要知道[2, sqrt(N)]范围内的所有素数，而这些素数恰好可能在进程 0 的区间里。所以并行化的第一个约束就是：每个进程都必须拥有完整的基素数表，基素数表的大小是sqrt(N)。当N = 10^10时，sqrt(N) = 10^5，这个表很小，复制到每个进程完全可行。真正的问题在于：进程 1 的区间里，哪些数是合数？它们的最小素因子可能小于N/2，也可能大于N/2。如果最小素因子在进程 0 的区间里，进程 0 必须把标记信息传给进程 1，否则进程 1 就会漏标。

2.2 按块划分 vs 按轮划分：两种主流策略的取舍

常见的 MPI 筛法划分策略有两种。第一种是按块划分：把[2, N]均匀切成P段，每个进程负责一段。每个进程先用串行筛法筛出[2, sqrt(N)]内的所有素数，然后对于自己区间内的每个基素数p，计算它在自己区间内的第一个倍数，从这个倍数开始标记。标记完成后，每个进程把自己区间内的素数收集起来，通过MPI_Allgatherv汇总。这种策略的优点是通信量小——只需要在最后汇总一次素数列表。缺点是负载不均衡：小区间内的素数密度高，标记操作多；大区间内素数密度低，但区间长度大，标记的总次数可能差不多，但缓存局部性差。

第二种是按轮划分：每一轮处理一个基素数p，所有进程协同标记p的所有倍数。进程r负责标记[r * chunk, (r+1) * chunk]范围内p的倍数。这种策略的优点是负载均衡极好，每一轮的标记量几乎相等。缺点是通信频繁：每一轮都需要广播当前的基素数p，如果sqrt(N) = 10^5，就有将近 10^4 个素数需要广播，通信开销会吃掉大部分并行收益。

我一般会选按块划分，因为它的通信模式简单，而且在实际测试中，当N大于10^9时，计算时间远大于通信时间。按轮划分更适合N较小但进程数很多的场景，比如N = 10^7、P = 64，这时候计算量太小，通信延迟反而成了瓶颈。

2.3 基素数表的获取：每个进程都要算一遍吗

基素数表是[2, sqrt(N)]内的所有素数。最直接的做法是每个进程独立跑一遍串行筛法，筛出sqrt(N)内的素数。当sqrt(N) = 10^5时，这个计算量微乎其微，完全没必要让进程 0 算完再广播。但如果N更大，比如N = 10^12，sqrt(N) = 10^6，串行筛10^6仍然很快，独立计算依然划算。真正需要广播基素数表的场景是：基素数表本身很大，或者生成基素数表的代价很高（比如需要试除到sqrt(sqrt(N))）。在埃拉托斯特尼筛法的语境下，基素数表永远很小，独立计算是最优解。

这里有一个容易翻车的点：基素数表的上界是sqrt(N)，不是N。如果你不小心把上界设成了N，每个进程都会试图分配一个大小为N的布尔数组，内存直接爆掉。血泪经验：在代码里把sqrt(N)单独算出来，存成一个变量，所有跟基素数表相关的循环都用这个变量做上界。

3. 用 MPI 实现按块划分筛法：从环境搭建到跑通第一个版本

3.1 环境准备与编译命令

假设你已经装好了 MPICH 或 OpenMPI，用mpicc编译。先写一个最小可运行的版本，验证通信和划分逻辑是否正确。创建文件mpi_sieve.c：

#include <mpi.h> #include <stdio.h> #include <stdlib.h> #include <math.h> #include <string.h> // 串行筛法，筛出 [2, limit] 内的素数，返回素数个数，primes 存储素数 int serial_sieve(long long limit, int *primes) { char *composite = (char *)calloc(limit + 1, sizeof(char)); int count = 0; for (long long i = 2; i <= limit; i++) { if (!composite[i]) { primes[count++] = (int)i; if (i * i <= limit) { for (long long j = i * i; j <= limit; j += i) { composite[j] = 1; } } } } free(composite); return count; } int main(int argc, char *argv[]) { MPI_Init(&argc, &argv); int rank, size; MPI_Comm_rank(MPI_COMM_WORLD, &rank); MPI_Comm_size(MPI_COMM_WORLD, &size); long long N = 1000000000LL; // 默认筛 10^9 if (argc > 1) N = atoll(argv[1]); long long sqrtN = (long long)sqrt((double)N); while (sqrtN * sqrtN > N) sqrtN--; while ((sqrtN + 1) * (sqrtN + 1) <= N) sqrtN++; // 每个进程独立计算基素数表 int *base_primes = (int *)malloc((sqrtN + 1) * sizeof(int)); int base_count = serial_sieve(sqrtN, base_primes); // 按块划分区间 [2, N] long long total = N - 1; // 2 到 N 共 N-1 个数 long long block = (total + size - 1) / size; long long start = 2 + rank * block; long long end = start + block - 1; if (end > N) end = N; if (start > N) { start = N + 1; end = N; } // 空区间 long long local_len = end - start + 1; char *local_composite = (char *)calloc(local_len > 0 ? local_len : 1, sizeof(char)); // 标记本地区间内的合数 for (int i = 0; i < base_count; i++) { long long p = base_primes[i]; if (p * p > N) break; // 计算 p 在本地区间内的第一个倍数 long long first = ((start + p - 1) / p) * p; if (first < p * p) first = p * p; for (long long j = first; j <= end; j += p) { local_composite[j - start] = 1; } } // 统计本地素数个数 long long local_count = 0; for (long long i = 0; i < local_len; i++) { if (!local_composite[i]) local_count++; } // 汇总所有进程的素数个数 long long total_count = 0; MPI_Reduce(&local_count, &total_count, 1, MPI_LONG_LONG, MPI_SUM, 0, MPI_COMM_WORLD); if (rank == 0) { printf("N = %lld, 素数总数 = %lld\n", N, total_count); } free(base_primes); free(local_composite); MPI_Finalize(); return 0; }

编译命令：

mpicc -O3 -o mpi_sieve mpi_sieve.c -lm

运行命令（4 个进程）：

mpirun -np 4 ./mpi_sieve 1000000000

这段代码的逻辑说明：每个进程先独立筛出[2, sqrt(N)]内的基素数，然后根据rank计算自己负责的区间[start, end]。对于每个基素数p，计算它在本地区间内的第一个倍数first，注意first不能小于p * p，因为小于p * p的倍数已经被更小的素因子标记过了。然后从first开始，每隔p标记一个合数。最后用MPI_Reduce汇总素数个数。

参数说明：N是筛法上界，通过命令行传入；sqrtN是基素数表上界，用整数平方根计算，注意浮点精度问题，需要做修正；block是每个进程负责的区间长度，用(total + size - 1) / size向上取整，保证所有数都被覆盖；local_composite是本地标记数组，下标i对应数字start + i。

3.2 通信模式的选择：Allgatherv 还是 Reduce

上面的代码只汇总了素数个数，没有汇总素数列表。如果你需要所有素数，比如做后续的密码学计算，就需要用MPI_Allgatherv把每个进程的素数收集起来。每个进程先把自己区间内的素数存到一个本地数组，然后调用MPI_Allgatherv。这里有一个坑：每个进程的素数个数不一样，所以需要先MPI_Allgather收集所有进程的素数个数，再MPI_Allgatherv收集素数本身。

// 收集本地素数 int *local_primes = (int *)malloc(local_count * sizeof(int)); long long idx = 0; for (long long i = 0; i < local_len; i++) { if (!local_composite[i]) local_primes[idx++] = (int)(start + i); } // 收集所有进程的素数个数 int *counts = (int *)malloc(size * sizeof(int)); MPI_Allgather(&local_count, 1, MPI_INT, counts, 1, MPI_INT, MPI_COMM_WORLD); // 计算位移 int *displs = (int *)malloc(size * sizeof(int)); displs[0] = 0; for (int i = 1; i < size; i++) { displs[i] = displs[i-1] + counts[i-1]; } // 收集所有素数 int *all_primes = (int *)malloc(total_count * sizeof(int)); MPI_Allgatherv(local_primes, local_count, MPI_INT, all_primes, counts, displs, MPI_INT, MPI_COMM_WORLD);

逻辑说明：MPI_Allgather收集每个进程的local_count，存到counts数组；displs是每个进程的素数在all_primes中的起始位置；MPI_Allgatherv把每个进程的local_primes按counts和displs拼接到all_primes。参数说明：local_count是本地素数个数，类型是long long，但MPI_Allgather要求int，所以需要强制转换；total_count是所有进程素数个数之和，用来分配all_primes。

注意：MPI_Allgatherv的counts和displs参数在非根进程上也会被使用，所以所有进程都必须正确计算这两个数组。不要只在 rank 0 上计算。

3.3 性能优化的第一步：减少标记次数

上面的代码有一个明显的性能问题：对于每个基素数p，都从p * p开始标记。但实际上，如果p * p小于start，那么first会被修正为start以上的第一个倍数。这个修正逻辑是对的，但标记的总次数并没有减少。真正能减少标记次数的优化是：只对奇数进行标记。因为除了 2 以外，所有偶数都是合数。所以可以把区间内的偶数直接跳过，只标记奇数。这样标记次数直接减半。

// 只标记奇数 for (int i = 1; i < base_count; i++) { // 跳过 base_primes[0] = 2 long long p = base_primes[i]; if (p * p > N) break; long long first = ((start + p - 1) / p) * p; if (first < p * p) first = p * p; if (first % 2 == 0) first += p; // 保证 first 是奇数 for (long long j = first; j <= end; j += 2 * p) { local_composite[j - start] = 1; } }

逻辑说明：跳过p = 2，因为偶数已经在初始化时被标记了。对于奇数p，first如果是偶数，就加上p变成奇数。然后步长变成2 * p，这样只标记奇数倍数。参数说明：base_primes[0]是 2，从i = 1开始循环；first % 2 == 0判断first是否为偶数，如果是，加上p变成奇数。

这个优化在N = 10^9时，标记次数从大约N * log(log(N))降到N/2 * log(log(N))，实测耗时减少约 40%。但要注意：初始化local_composite的时候，需要把偶数对应的位置直接标记为 1，否则偶数会被误判为素数。

4. 避坑与排查：MPI 筛法最容易翻车的五个地方

4.1 现象：筛出来的素数个数比串行版本少

原因：区间划分的边界计算错误。最常见的是block计算用了整数除法，导致最后一个进程的区间没有覆盖到N。比如total = N - 1，size = 4，block = (N-1+3)/4，如果N-1不能被 4 整除，最后一个进程的end会小于N。解决：在计算end之后，强制if (end > N) end = N;，并且确保start不超过N。另外，如果start > N，说明这个进程没有任务，local_len设为 0，不要分配负数长度的数组。

4.2 现象：程序在MPI_Allgatherv处挂死

原因：counts和displs数组在非根进程上没有正确初始化。MPI_Allgatherv要求所有进程都提供完整的counts和displs，但很多人在写代码时只在 rank 0 上计算了这两个数组。解决：在所有进程上都执行MPI_Allgather和displs计算，不要用if (rank == 0)包起来。另外，counts的类型是int，但local_count是long long，如果local_count超过INT_MAX，会溢出。对于N = 10^10，素数个数大约4.5 * 10^8，不会溢出，但N = 10^12时就要小心了。

4.3 现象：标记数组越界，程序崩溃

原因：local_composite的下标计算错误。local_composite[j - start]中，j的范围是[first, end]，start是区间起点，所以j - start的范围是[0, end - start]。但如果first小于start，j - start会变成负数。解决：在计算first之后，强制if (first < start) first = start;，并且确保first是p的倍数。另外，local_len的计算是end - start + 1，如果start > end，local_len为负数，calloc会失败。解决：在start > N时，把local_len设为 0，local_composite分配 1 个字节。

4.4 现象：多进程运行时结果不稳定，有时对有时错

原因：MPI_Reduce的MPI_LONG_LONG类型和local_count的类型不匹配。local_count是long long，但MPI_Reduce的count参数是1，datatype是MPI_LONG_LONG，这本身没问题。但如果local_count被声明为int，而MPI_Reduce用了MPI_LONG_LONG，就会读错内存。解决：统一用long long声明local_count和total_count，并且MPI_Reduce的datatype用MPI_LONG_LONG。另外，MPI_Reduce的op用MPI_SUM，不要用MPI_MAX或MPI_MIN。

4.5 现象：性能没有随进程数增加而提升，反而变慢

原因：通信开销超过了计算收益。当N较小（比如N = 10^7）而进程数较多（比如P = 32）时，每个进程的计算量很小，但MPI_Allgatherv的通信量很大。解决：对于小N，减少进程数，或者改用按轮划分策略，减少通信次数。另外，检查是否开启了编译优化，-O3是必须的。如果用了-O0，串行代码本身就很慢，并行收益会被掩盖。还有一个常见原因是：基素数表在每个进程上独立计算，如果sqrt(N)很大，这个计算量不可忽略。解决：让 rank 0 计算基素数表，然后用MPI_Bcast广播给其他进程。但如前所述，对于N = 10^9，sqrt(N) = 31623，独立计算完全没问题。

5. 进阶技巧：用 MPI-3 共享内存和非阻塞通信再压榨一轮

5.1 用MPI_Win_allocate_shared减少内存复制

MPI-3 引入了共享内存窗口，允许同一节点上的进程直接访问彼此的内存。对于筛法来说，基素数表是所有进程都要读的，如果每个进程都存一份，内存浪费不说，缓存命中率也低。用MPI_Win_allocate_shared可以让所有进程共享同一份基素数表。

MPI_Win win; int *shared_base; MPI_Win_allocate_shared(sqrtN * sizeof(int), sizeof(int), MPI_INFO_NULL, MPI_COMM_WORLD, &shared_base, &win); // rank 0 写入基素数表 if (rank == 0) { int *base_primes = (int *)malloc(sqrtN * sizeof(int)); int base_count = serial_sieve(sqrtN, base_primes); memcpy(shared_base, base_primes, base_count * sizeof(int)); free(base_primes); } MPI_Win_fence(0, win); // 同步 // 所有进程现在可以读 shared_base

逻辑说明：MPI_Win_allocate_shared在每个进程上分配一段共享内存，shared_base是指向本地共享内存的指针。rank 0 把基素数表写入自己的共享内存，MPI_Win_fence同步后，其他进程可以通过MPI_Win_shared_query获取 rank 0 的共享内存指针，直接读取。参数说明：sqrtN * sizeof(int)是共享内存大小，sizeof(int)是每个元素的字节数，MPI_INFO_NULL表示不使用额外信息。注意：MPI_Win_fence是集体操作，所有进程都必须调用。

这个优化在单节点多进程场景下效果明显，因为基素数表只存了一份，L3 缓存命中率提高。但在多节点场景下，跨节点的共享内存不可用，还是需要MPI_Bcast。

5.2 用非阻塞通信重叠计算和通信

按块划分的筛法在最后汇总素数列表时，MPI_Allgatherv是阻塞的。如果素数列表很大，通信时间会很长。可以用MPI_Iallgatherv非阻塞版本，在通信的同时做一些收尾计算，比如统计素数个数、计算校验和。

MPI_Request req; MPI_Iallgatherv(local_primes, local_count, MPI_INT, all_primes, counts, displs, MPI_INT, MPI_COMM_WORLD, &req); // 在等待通信完成的同时，计算本地素数的校验和 long long local_checksum = 0; for (long long i = 0; i < local_count; i++) { local_checksum += local_primes[i]; } MPI_Wait(&req, MPI_STATUS_IGNORE);

逻辑说明：MPI_Iallgatherv立即返回一个MPI_Request，通信在后台进行。在MPI_Wait之前，可以执行不依赖all_primes的计算。参数说明：req是请求句柄，MPI_Wait等待通信完成。注意：MPI_Iallgatherv的缓冲区all_primes在MPI_Wait返回之前不能被修改，否则会读到脏数据。

5.3 验证方法：用已知素数个数做交叉验证

筛法结果对不对，最直接的验证方法是查素数计数函数π(N)的已知值。比如π(10^9) = 50847534，π(10^10) = 455052511。跑完程序后，把输出的素数个数和这些已知值对比。如果一致，说明筛法逻辑正确。如果不一致，先检查区间划分是否覆盖了[2, N]的所有数，再检查基素数表是否完整，最后检查标记逻辑是否有遗漏。

另一个验证方法是：用串行版本跑一个小N（比如N = 10^6），把素数列表存到文件，然后用 MPI 版本跑同样的N，对比两个文件的素数列表是否完全一致。这个方法能定位到具体是哪个素数漏了或者多了。

5.4 一个具体技巧：用位数组压缩内存

local_composite用的是char数组，每个元素占 1 字节。如果改成位数组，每个元素只占 1 位，内存占用直接降到 1/8。对于N = 10^10，char数组需要 10 GB，位数组只需要 1.25 GB。实现方法是用uint64_t数组，每个uint64_t存 64 个数的标记。

uint64_t *bits = (uint64_t *)calloc((local_len + 63) / 64, sizeof(uint64_t)); // 设置第 i 位 void set_bit(uint64_t *bits, long long i) { bits[i / 64] |= (1ULL << (i % 64)); } // 读取第 i 位 int get_bit(uint64_t *bits, long long i) { return (bits[i / 64] >> (i % 64)) & 1; }

逻辑说明：bits数组的每个元素是uint64_t，能存 64 个标记。set_bit把第i位设为 1，get_bit读取第i位。参数说明：(local_len + 63) / 64是向上取整计算需要的uint64_t个数。注意：位数组的标记和读取比char数组慢，因为需要位运算。实测在N = 10^9时，位数组版本比char数组版本慢约 15%，但内存占用降到 1/8。如果内存是瓶颈，位数组值得用；如果 CPU 是瓶颈，还是用char数组。

我自己的习惯是：先用char数组跑通逻辑，确认结果正确后，再根据内存和性能的权衡决定是否换成位数组。不要一上来就搞位数组，调试难度会大很多。希望帮到你。

本文还有配套的精品资源，点击获取