c - OpenMP C 正态随机发生器

Question

我正在为并行计算进行 C 作业，其中我必须使用 Box-Muller 变换通过高效的安全正态随机生成器来实现某种蒙特卡罗模拟。我生成 2 个均匀随机数 X 和 Y 的 vector ，条件是 X 在 (0,1] 中，Y 在 [0,1] 中。但我不确定我从半开区间 (0,1] 中采样均匀随机数的方法是否正确。有人遇到过类似的事情吗？我正在使用以下代码:

double* StandardNormalRandom(long int N){
    double *X = NULL, *Y = NULL, *U = NULL;

    X = vUniformRandom_0(N / 2);
    Y = vUniformRandom(N / 2);
    #pragma omp parallel for
    for (i = 0; i<N/2; i++){
            U[2*i] = sqrt(-2 * log(X[i]))*sin(Y[i] * 2 * pi);
            U[2*i + 1] = sqrt(-2 * log(X[i]))*cos(Y[i] * 2 * pi);
            }
    return U;
    }

    double* NormalRandom(long int N, double mu, double sigma2)
    {
            double *U = NULL, stdev = sqrt(sigma2);
            U = StandardNormalRandom(N);
            #pragma omp parallel for
            for (int i = 0; i < N; i++) U[i] = mu + stdev*U[i];
            return U;
    }

这是我的 UniformRandom 函数的一部分，也是并行实现的:

#pragma omp parallel for  firstprivate(i)
for (long int j = 0; j < N;j++)
{
    if (i == 0){
        int tn = omp_get_thread_num();
        I[tn] = S[tn];
        i++;
    }
    else
    {
        I[j] = (a*I[j - 1] + c) % m;
    }

}

}
#pragma omp parallel for 
for (long int j = 0; j < N; j++)
    U[j] = (double)I[j] / (m+1.0);

Answer 1

在 StandardNormalRandom 函数中，我假设指针 U 已分配给大小 N，在这种情况下，该函数看起来对我来说很好。以及函数 NormalRandom。

但是对于函数 UniformRandom (缺少一些部分，所以我必须假设一些东西)，如果以下行 I[j] = (a*I[ j - 1] + c) % m + 1; 是带有 omp parallel for 的循环体，那么您将会遇到一些问题。由于无法知道线程的执行顺序，因此当前线程(固定值为j)无法依赖I[j - 1]<的值 因为该值可以随时修改(I 应默认共享)。

希望对你有帮助!