如何用OpenMP多线程改写CPILL程序，加速π值计算？

一、OpenMP性能倒退的原因 #

• 直接对全局累加变量进行多线程读写会导致缓存乒乓（Cache Ping-Pong）：多个核心的缓存不断同步累加变量的状态，开销远大于多线程带来的并行收益。
• MPFR高精度变量的操作本身开销较大，全局共享时的锁竞争或内存同步成本被进一步放大。

二、改写OpenMP优化代码的思路 #

核心是分块独立计算+最后合并，每个线程维护自己的MPFR累加变量，彻底避免跨线程竞争：

#include <mpfr.h>
#include <omp.h>
#include <stdio.h>
#include <stdlib.h>

int main(int argc, char *argv[]) {
    if (argc != 2) {
        fprintf(stderr, "Usage: %s <iterations>\n", argv[0]);
        return 1;
    }
    unsigned long long total_iter = strtoull(argv[1], NULL, 10);
    int num_threads = omp_get_max_threads();
    printf("Using %d threads\n", num_threads);

    // 为每个线程分配独立的MPFR累加变量
    mpfr_t *thread_sums = malloc(num_threads * sizeof(mpfr_t));
    mpfr_t pi;
    mpfr_init2(pi, 1000); // 设置精度为1000位（可根据需求调整）
    for (int i = 0; i < num_threads; i++) {
        mpfr_init2(thread_sums[i], 1000);
        mpfr_set_ui(thread_sums[i], 0, MPFR_RNDN);
    }

    #pragma omp parallel
    {
        int tid = omp_get_thread_num();
        // 分配每个线程的迭代区间，最后一个线程处理剩余部分
        unsigned long long start = tid * (total_iter / num_threads);
        unsigned long long end = (tid == num_threads - 1) ? total_iter : (tid + 1) * (total_iter / num_threads);
        
        mpfr_t term, denom;
        mpfr_init2(term, 1000);
        mpfr_init2(denom, 1000);

        for (unsigned long long i = start; i < end; i++) {
            mpfr_set_ui(denom, 2 * i + 1, MPFR_RNDN);
            mpfr_ui_div(term, 1, denom, MPFR_RNDN);
            // 根据迭代项奇偶性加减
            if (i % 2 == 0) {
                mpfr_add(thread_sums[tid], thread_sums[tid], term, MPFR_RNDN);
            } else {
                mpfr_sub(thread_sums[tid], thread_sums[tid], term, MPFR_RNDN);
            }
        }

        mpfr_clear(term);
        mpfr_clear(denom);
    }

    // 合并所有线程的计算结果
    mpfr_set_ui(pi, 0, MPFR_RNDN);
    for (int i = 0; i < num_threads; i++) {
        mpfr_add(pi, pi, thread_sums[i], MPFR_RNDN);
        mpfr_clear(thread_sums[i]);
    }
    mpfr_mul_ui(pi, pi, 4, MPFR_RNDN); // Leibniz公式结果乘以4得到π

    mpfr_printf("π = %.50Rf\n", pi);

    mpfr_clear(pi);
    free(thread_sums);
    return 0;
}

编译命令（需确保已安装MPFR和GMP库）：

gcc -O3 -fopenmp cpill.c -o cpill -lmpfr -lgmp

代码优化点说明 #

• 每个线程拥有独立的累加变量，彻底消除跨线程内存竞争。
• 分块时确保最后一个线程处理剩余迭代，避免计算遗漏。
• 提前初始化MPFR变量并指定精度，减少运行时动态内存开销。

三、其他可行的加速方案 #

1. 替换为收敛速度更快的公式 #

Leibniz公式收敛极慢（每约4次迭代仅提升1位十进制精度），1e15次迭代的性价比极低。推荐使用Chudnovsky公式，每计算一项可提升约14位精度，计算千位精度仅需几十项，效率远超Leibniz公式。

2. SIMD指令优化 #

利用AMD Ryzen 5 3600支持的AVX2指令集，对低精度中间项进行向量并行计算。由于MPFR直接SIMD支持有限，可先用原生浮点指令计算块内近似和，再用MPFR进行高精度修正，平衡速度与精度。

3. 分布式计算 #

若必须使用Leibniz公式完成1e15次迭代，可将任务拆分到多台机器，每台机器计算一部分迭代的和，最后合并结果。你的64GB内存足够存储大量中间结果，适合分布式分块计算。

4. 内存与缓存优化 #

• 使用posix_memalign分配对齐内存，将MPFR变量固定在CPU的L2/L3缓存中，减少内存访问延迟。
• 降低MPFR变量精度到实际需求水平：若不需要上万位精度，降低精度可大幅提升计算速度。

5. 专用高精度计算库替代 #

考虑使用GMP整数运算结合分数累加（Leibniz本质是分数求和），或参考y-cruncher这类π计算专用工具的极致优化思路，这类工具针对大规模π计算做了多层并行、缓存优化等深度调优。

内容的提问来源于stack exchange，提问作者Sylwester Bogusiak