OpenMP并行蒙特卡洛积分异常：无printf结果错，有printf结果正确

问题描述 #

尝试用蒙特卡洛方法计算区间[a,b]内的定积分，单线程运行正常，但添加#pragma omp parallel for并行指令后出现异常：注释掉循环内的printf语句时，计算结果错误（例如a=0、b=1时，正确结果应为1/3）；保留printf语句时结果却正确。运行环境为Windows 10 + Visual Studio，代码如下：

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

float FloatRandomizer(float a, float b)
{
    float randomNumber = (a + (float)rand() / (float)(RAND_MAX / b));

    if (randomNumber > b)
        randomNumber = randomNumber - a;

    return randomNumber;
}

float F(float x)
{
    return x * x;
}


int main()
{
    int n = 1000000;

    float a;
    printf("a = "); scanf_s("%f", &a);

    float b;
    printf("b = "); scanf_s("%f", &b);

    float temp = (b - a) / n;

    float sum = 0;
    int i;

#pragma omp parallel for shared(sum)
    for (i = 0; i < n; i++)
    {
        float randomNumber = FloatRandomizer(a, b);
        sum = sum + F(randomNumber);
        //printf("threadNum = %d\nsum = %f\n", omp_get_thread_num(), sum);
    }

    float result = temp * sum;
    printf("result = %f", result);
}

问题原因 #

• 竞态条件导致sum计算错误：多个线程同时对共享变量sum执行sum = sum + F(randomNumber)操作，该操作分为"读取sum→计算新值→写回sum"三步，并非原子操作。若线程A读取sum后还未写回，线程B就读取了同一个sum值，两者计算后写回会丢失其中一个线程的贡献，最终sum值偏小。
• printf的隐式同步掩盖问题：printf是耗时操作且带有隐式同步，会减慢线程执行速度，大幅降低多线程同时操作sum的概率，相当于无意中给sum的更新加了同步，因此结果看似正确，但这只是巧合，并非正确解决方案。
• rand()的线程不安全问题：rand()函数不是线程安全的，多线程同时调用会导致随机数序列混乱，进一步影响积分结果的准确性。

解决方法 #

1. 用Reduction消除竞态条件
   最高效的方式是使用OpenMP的reduction子句，让每个线程拥有独立的sum私有副本，循环结束后自动合并所有副本的值到全局sum，避免线程间的竞争：

   float sum = 0;
   int i;
   
   #pragma omp parallel for reduction(+:sum)
   for (i = 0; i < n; i++)
   {
       float randomNumber = FloatRandomizer(a, b);
       sum += F(randomNumber);
   }
   

   也可以用原子操作临时保护sum的更新，但效率低于reduction：

   #pragma omp parallel for shared(sum)
   for (i = 0; i < n; i++)
   {
       float randomNumber = FloatRandomizer(a, b);
       #pragma omp atomic
       sum += F(randomNumber);
   }
   

2. 替换为线程安全的随机数生成
   在Windows环境下，用rand_s()替代rand()，保证多线程下随机数生成的安全性：

   #include <windows.h> // 需要包含头文件
   
   float FloatRandomizer(float a, float b)
   {
       unsigned int randomValue;
       rand_s(&randomValue);
       // 直接生成[a,b]区间的随机数，无需额外判断
       float randomNumber = a + (float)randomValue / UINT_MAX * (b - a);
       return randomNumber;
   }
   

3. 简化积分公式（可选）
   原代码的积分公式计算逻辑正确，可简化写法让逻辑更清晰：

   float result = (b - a) * sum / n;
   

内容的提问来源于stack exchange，提问作者Good guy from world