为何OpenMP并行蒙特卡洛PI估算结果与预期不符？

问题描述 #

从YouTube视频及橡树岭国家实验室官网获取的OpenMP并行蒙特卡洛PI估算代码，在VS Studio 22、Windows 11环境（16核i9-12900kf、32GB内存）运行时，得到的PI值仅为预期值的约一半。将第二个代码的随机函数改为srand()和rand()后，问题依旧。

原始代码 #

代码1 #

int main() {
    int num; // number of iterations
    printf("Enter number of iterations you want the loop to run for:   ");
    scanf_s("%d", &num);

    double x, y, z, pi;
    long long int i;
    int count = 0;
    int num_thread;

    printf("Enter number of threads you want to run to parallelize the process:\t");
    scanf_s("%d", &num_thread);
    printf("\n");

#pragma omp parallel firstprivate(x,y,z,i) shared(count) num_threads(num_thread)
    {
        srand((int)time(NULL) ^ omp_get_thread_num());
        for (i = 0; i < num; i++) {
            x = (double)rand() / (double)RAND_MAX;
            y = (double)rand() / (double)RAND_MAX;

            z = pow(((x * x) + (y * y)), .5);

            if (z <= 1) {
                count++;
            }
        }
    } // END PRAGMA
    pi = ((double)count / (double)(num * num_thread)) * 4;
    printf("The value of pi obtained is %f\n", pi);

    return 0;
}

代码2（已修改为srand()/rand()） #

#include <stdio.h>
#include <stdlib.h>
#include <omp.h>
#include <math.h>
 
int main(int argc, char* argv[])
{
    int niter = 1000000;            //number of iterations per FOR loop
    double x,y;                     //x,y value for the random coordinate
    int i;                          //loop counter
    int count=0;                //Count holds all the number of how many good coordinates
    double z;                       //Used to check if x^2+y^2<=1
    double pi;                      //holds approx value of pi
    int numthreads = 16;
 
    #pragma omp parallel firstprivate(x, y, z, i) shared(count) num_threads(numthreads)
    {
        srand((int)time(NULL) ^ omp_get_thread_num());    //Give random() a seed value
        for (i=0; i<niter; ++i)              //main loop
        {
            x = (double)rand()/RAND_MAX;      //gets a random x coordinate
            y = (double)rand()/RAND_MAX;      //gets a random y coordinate
            z = sqrt((x*x)+(y*y));          //Checks to see if number is inside unit circle
            if (z<=1)
            {
                ++count;            //if it is, consider it a valid random point
            }
        }
        //print the value of each thread/rank
    }
    pi = ((double)count/(double)(niter*numthreads))*4.0;
    printf("Pi: %f\n", pi);
 
    return 0;
}

问题根源 #

绝非机器硬件问题，核心原因是两个关键bug：

1. 共享变量count的线程竞争：count++是包含读取-修改-写入的非原子操作，多个线程同时操作时会出现数据丢失。例如两个线程同时读取count=100，各自加1后写回，最终count仅变为101而非预期的102，导致总计数远低于实际有效点数，PI值被低估。
2. 随机数种子冲突：srand((int)time(NULL) ^ omp_get_thread_num())依赖秒级时间戳，若线程启动间隔小于1秒，不同线程的种子会重复，导致随机数序列完全一致，采样点分布重复，进一步加剧结果偏差。

修复方案 #

方案1：原子操作保护共享计数 #

将count++;替换为原子操作，确保每次自增的完整性：

#pragma omp atomic
count++;

方案2：私有计数器+全局合并（更高效） #

每个线程维护私有计数器，并行结束后再合并结果，减少共享变量竞争：

#pragma omp parallel firstprivate(x,y,z,i) num_threads(num_thread)
{
    // 用高精度时间戳生成唯一种子
    unsigned int seed = (unsigned int)(omp_get_wtime() * 1000000) ^ omp_get_thread_num();
    srand(seed);
    
    int local_count = 0; // 线程私有计数器
    for (i = 0; i < num; i++) {
        x = (double)rand() / (double)RAND_MAX;
        y = (double)rand() / (double)RAND_MAX;
        z = sqrt(x*x + y*y); // 替代pow提升效率
        if (z <= 1) {
            local_count++;
        }
    }
    // 原子操作合并私有计数到全局
    #pragma omp atomic
    count += local_count;
}

修复后的完整代码（以代码1为例） #

#include <stdio.h>
#include <stdlib.h>
#include <omp.h>
#include <math.h>
#include <time.h>

int main() {
    int num; // 迭代次数
    printf("请输入循环迭代次数:   ");
    scanf_s("%d", &num);

    double x, y, z, pi;
    long long int i;
    int count = 0;
    int num_thread;

    printf("请输入并行化使用的线程数:\t");
    scanf_s("%d", &num_thread);
    printf("\n");

#pragma omp parallel firstprivate(x,y,z,i) num_threads(num_thread)
    {
        // 生成唯一的线程随机种子
        unsigned int seed = (unsigned int)(omp_get_wtime() * 1000000) ^ omp_get_thread_num();
        srand(seed);
        
        int local_count = 0;
        for (i = 0; i < num; i++) {
            x = (double)rand() / (double)RAND_MAX;
            y = (double)rand() / (double)RAND_MAX;
            z = sqrt(x*x + y*y);
            
            if (z <= 1) {
                local_count++;
            }
        }
        // 原子操作合并结果
        #pragma omp atomic
        count += local_count;
    }

    pi = ((double)count / (double)(num * num_thread)) * 4;
    printf("估算得到的PI值为: %f\n", pi);

    return 0;
}

关键修复点说明 #

• 用omp_get_wtime()的高精度时间戳替代time(NULL)，确保每个线程的随机数种子唯一，避免序列重复。
• 线程私有计数器减少共享变量的竞争次数，提升并行效率。
• 原子操作保证全局计数合并的线程安全，避免数据丢失。
• 用sqrt()替代pow()计算模长，性能更优且结果一致。

内容的提问来源于stack exchange，提问作者Adam G