OpenMP并行区域for循环被忽略问题求助（蒙特卡洛求Pi）

Let's walk through the main issues in your code, starting with why the loop is being ignored when using multiple threads:

1. Integer Division Kills Your Loop Iterations #

The biggest problem is this line:

for(int i=0;i<(num_points*(1/nthreads));i++)

Since nthreads is an integer, 1/nthreads uses integer division. When nthreads >= 2, 1/nthreads evaluates to 0 (because integer division truncates towards zero). Multiplying that by num_points gives you 0 iterations—so the loop never runs.

To fix this, use floating-point division (cast one operand to double) or better yet, use integer arithmetic to split the workload evenly (and handle any leftover points):

int points_per_thread = num_points / nthreads;
int remainder = num_points % nthreads;
// Give the first 'remainder' threads one extra point to avoid losing total count
int my_points = points_per_thread + (omp_get_thread_num() < remainder ? 1 : 0);

2. Thread-Unsafe Random Number Generation #

Your current random number setup has two critical flaws:

• Duplicate seeds: Multiple threads will likely call time(0) at almost the same time, getting the same timestamp. Even multiplying by tid_set_clock might not save you here—if the timestamp is identical, seeds could still collide.
• Non-thread-safe rand(): The standard rand() function isn't thread-safe. When multiple threads call it simultaneously, you'll get race conditions and corrupted random numbers.

Fix this by using a thread-safe random number generator like rand_r() (POSIX-compliant) with a thread-local seed:

unsigned int seed = (unsigned int)(time(NULL) ^ omp_get_thread_num()); // Unique seed per thread
// Inside the loop:
double coord[0] = 2.0 * ((double)rand_r(&seed) / RAND_MAX);
double coord[1] = 2.0 * ((double)rand_r(&seed) / RAND_MAX);

3. points_inside Isn't Being Aggregated #

Right now, points_inside is marked as private, meaning each thread has its own copy. But you never combine these values into a global count—so even if the loop ran, you'd lose all the thread-specific results.

Use OpenMP's reduction clause to automatically sum up each thread's count:

long long total_inside = 0;
#pragma omp parallel num_threads(threads) reduction(+:total_inside)
{
    int my_inside = 0;
    // ... loop logic ...
    if(dist_cent <= 1.0) { // Side note: Your original logic was reversed! Points inside have dist <=1
        my_inside++;
    }
    total_inside += my_inside;
}

(Also, your original code had dist_cent >=1 then points_inside--—that's backwards for counting points inside the circle. Fixed that above.)

Full Corrected Code Example #

Here's a cleaned-up version of your code with all fixes applied:

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

int main() {
    const int num_points = 1000000;
    const int threads = 4;
    long long total_inside = 0;

    #pragma omp parallel num_threads(threads) reduction(+:total_inside)
    {
        int tid = omp_get_thread_num();
        int nthreads = omp_get_num_threads();
        
        // Assign workload evenly, handle remainder
        int points_per_thread = num_points / nthreads;
        int remainder = num_points % nthreads;
        int my_points = points_per_thread + (tid < remainder ? 1 : 0);

        // Thread-safe random seed
        unsigned int seed = (unsigned int)(time(NULL) ^ tid);
        int my_inside = 0;

        for(int i = 0; i < my_points; i++) {
            double x = 2.0 * ((double)rand_r(&seed) / RAND_MAX);
            double y = 2.0 * ((double)rand_r(&seed) / RAND_MAX);
            double dist = sqrt(pow(x - 1.0, 2) + pow(y - 1.0, 2));
            
            if(dist <= 1.0) {
                my_inside++;
            }
        }

        total_inside += my_inside;
    }

    double pi_estimate = 4.0 * ((double)total_inside / num_points);
    printf("Estimated Pi: %.6f\n", pi_estimate);
    return 0;
}

Compilation Note #

Don't forget to compile with OpenMP enabled:

gcc -fopenmp pi.c -o pi -lm

内容的提问来源于stack exchange，提问作者Armando Canales Lima