Estimating Pi: A Monte Carlo Approach Enhanced by OpenMP Parallelism

Image credit: domin_domin via iStock

What if I told you that the seemingly random nature of points can be harnessed to calculate one of mathematics’ most famous constants — π? In this blog, we will explore how the Monte Carlo method can be used to estimate the value of π, an essential constant in mathematics and physics. We will also examine how OpenMP parallelism enhances the efficiency of the Monte Carlo approach.

Understanding the Monte Carlo Method

Monte Carlo methods are a broad class of computational algorithms that rely on repeated random sampling to obtain numerical results. They are widely used in various fields, including statistics, finance, physics, and engineering, to solve problems that may be deterministic in nature but are difficult to solve analytically.

One interesting application of the Monte Carlo method is the estimation of π. By randomly generating points within a unit square and determining how many of those points fall within a quarter circle inscribed in that square, we can derive an estimate for π.

Algorithm for Estimating π

1. Draw a unit square and inscribe a quarter circle within it.
2. Randomly generate points within the square.
3. Count the number of points that fall within the quarter circle by checking whether x²+y²≤1.
4. Calculate the ratio of points inside the circle to the total points, which approximates π/4.
5. Estimate π by multiplying the ratio by 4.

Monte Carlo method applied to approximating the value of π

Parallel Implementation of Monte Carlo Method

By distributing the workload across multiple threads, the parallel implementation of the Monte Carlo method using OpenMP significantly reduces computation time, especially when handling large datasets. This approach allows for scaling the number of points generated, achieving a more accurate approximation of π. Below is the implementation of the parallel Monte Carlo method for estimating π using OpenMP.

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

int main() {
    long int n;
    long int count = 0;
    long double x, y, z;
    double start, end, wall_clock_time;

    printf("\nn = ");
    scanf("%ld", &n);

    omp_set_num_threads(8); // Set the number of threads

    // True random numbers from random.org (based on atmospheric noise)
    unsigned int true_random_numbers[8] = {74, 51, 54, 89, 54, 15, 23, 31};

    start = omp_get_wtime(); // Record the start time

    #pragma omp parallel
    {
        // Seed each thread once
        unsigned int seed = true_random_numbers[omp_get_thread_num()]; 
        long int local_count = 0;

        #pragma omp for private(x, y, z)
        for(long int i = 0; i < n; ++i) {
            x = (long double)rand_r(&seed) / RAND_MAX; // Generate random x in [0, 1]
            y = (long double)rand_r(&seed) / RAND_MAX; // Generate random y in [0, 1]
            z = x * x + y * y; // Check if inside the unit circle

            if(z <= 1) {
                local_count++; // Count points inside the circle
            }
        }

        #pragma omp atomic
        count += local_count;
    }

    end = omp_get_wtime(); // Record the end time

    long double pi = (long double)count / n * 4; // Estimate of Pi
    printf("\nApproximate PI = %.9Lf\n", pi);

    wall_clock_time = end - start;
    printf("Elapsed Wall Clock Time: %f seconds\n", wall_clock_time);

    return 0;
}

The provided C program estimates the value of π using the Monte Carlo method, enhanced with parallel processing via OpenMP. It uses true random numbers sourced from random.org, generated using atmospheric noise to ensure high-quality randomness. Each thread is seeded once using the true_random_numbers[] array, allowing for independent and thread-safe random number generation with rand_r().

The workload is effectively distributed among threads using the#pragma omp for directive, enabling each thread to execute the Monte Carlo algorithm on its assigned set of iterations. Each thread maintains a private copy of the local count of points that fall within the quarter circle. After all threads complete their assigned tasks, the local counts are combined using an atomic operation to prevent race conditions. Finally, the estimated value of π is calculated by multiplying the ratio of points inside the circle by 4.

Experimental Results: Estimating π with 80 Billion Points

In the experiment conducted to estimate the value of π using a substantial sample size of 80 billion points , the Monte Carlo method was executed with n=80,000,000,000 (80 billion), yielding an approximation of π=3.141592614, accurate up to 7 decimal places , and took an elapsed wall clock time of approximately 235.32 seconds ( about 3.92 minutes ) to approximate π.

For comparison, the serial version of the algorithm was run with n=8,000,000,000 ( 8 billion ), resulting in an approximation of π=3.14160, accurate up to 3 decimal places, and took an elapsed time of 147.69 seconds ( about 2.46 minutes ).

This showcases a significant efficiency gain with the parallel implementation; the serial version would require approximately 24 to 25 minutes to approximate π for a sample size of 80 billion.

Experiment Setup:

• Device: Apple MacBook Air
• Chip: M2
• RAM: 16GB
• CPU: 8-core (4 efficiency cores, 4 performance cores)
• Operating System: macOS Sonoma Version 14.6.1
• Compiler: GCC 14 with OpenMP support

Convergence of the Monte Carlo Estimation of Pi

The convergence graph shows how the sample size affects the estimated value of π using the Monte Carlo method. As the sample size increases, the estimated value of π stabilizes around the true value of π, which is indicated by a red dashed line. For instance, with a sample size of 10, the estimate is 3.6, while at 1,000,000 samples, it improves to about 3.139. This trend highlights that larger sample sizes yield more accurate estimates, confirming the effectiveness of the Monte Carlo method for approximating the value of Pi.

Absolute Error vs. Sample Size

The absolute error graph illustrates how the sample size impacts the accuracy of the π estimation using the Monte Carlo method. As the sample size increases, the absolute error decreases significantly, indicating that larger sample sizes lead to more precise approximations. For example, with a sample size of 10, the absolute error is 0.458407, whereas, at a sample size of 1,000,000, the error drops to just 0.002537. This trend demonstrates the reliability of the Monte Carlo method in achieving high accuracy with sufficient sample sizes, confirming its effectiveness for approximating the value of π.

Conclusion

In this study, we applied the Monte Carlo method to estimate π using a sample size of 80 billion points, yielding an approximation of 3.141592614, accurate up to 7 decimal places. The results demonstrated the effectiveness of parallel computing in reducing computation time compared to serial implementation. Additionally, the analysis of convergence and absolute error confirmed that larger sample sizes improve accuracy, demonstrating the robustness of the Monte Carlo method for π approximation.