Question

Develop a program using Threads in C/C++ to estimate the value of PI using the Monte Carlo method

use:

C/C++ #include srand((unsigned)(myid)); x = ((double)rand()) / ((double)RAND_MAX); y = ((double)rand()) / ((double)RAND_MAX);

Your program will allow the user to specify the number of threads (range 1 to 10) and the total number of data points (range 10 to 1,000,000) used for the Monte Carlo simulation on the command line. Note, DO NOT assume the number of data points is always evenly divisible by the number of threads.

Answer 1

There are many libraries for using threads in a C/C++ program. One of the most popular ones is pthreads. The following functions are used in pthreads:

pthread_create() :- Creates a thread, and calls a function to be executed by the thread. The argument of function is provided after the function to be executed. In the case of below program, the function's name is func. In the given program, pthread_create is called with 4 arguments. The first is the address of thread. The second argument is optional (so NULL is passed). The 3rd argument is the name of the function.The 4th argument is the parameter that is to be passed to the function. In this case, the argument is the id of the thread. The argument is passed as void pointer of long[Convention of pthreads].It returns 0 if thread was succesfully created. Otherwise it returns the error code.

pthread_exit() :- This is called to exit the pthread and its argument can be used in the pthread_join to store the output of the thread.

pthread_join() :- Terminates the thread.

Now about Monte Carlo algorithm. Monte Carlo algorithm finds the value of pi by finding points in a unit square(square of side length 1 inside the first quadrant) . A circle is drawn with radius 1. Then the point inside the square lies inside the circle if x^2 + y^2 < 1. Also 0<x<1 and 0<y<1. Area of the circle lying in first quadrant is the ratio of total number of points inside the circle to the ratio of total number of data points. Area of the whole circle = 4* Area of quadrant of circle. Also area of circle = pi*r*r= pi*1*1 = pi.

Thus pi = 4*(number of points inside the circle)/(total number of data points inside the square). Accuracy increases as the number of data points increase.

Since it is possible that the number of data points is not divisible by number of threads, it is required to efficiently distribute the tasks to threads. This is a design choice. Here assume that there are 3 threads and 20 data points. Then all the 3 threads will create 20/3 = 6 data points/thread. In addition, the first thread(thread number 0) executes the remaining points as well (20%3 = 2) .

The following program does the required tasks:(The program is in C language).

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

//Global variables
int num_threads, data_pts;

//Function executed by each thread
void *func(void * threadID) {
    long thid = (long)threadID;
    int tid = (int)thid;
    
    //Use malloc so that the value does not get destroyed on function return
    double * circ_count = (double *) malloc(sizeof(double));
    *circ_count = 0;
    
    int iter = data_pts/num_threads;
    //Make the thread do the required no of computations
    
    for(int c = 0;c<iter;c++) {
        double x = ((double)rand()) / ((double)RAND_MAX);
        double y = ((double)rand()) / ((double)RAND_MAX);
        double r = sqrt(x*x+y*y);
        if(r<1) *circ_count += 1;
    }
    
    //Make the first thread do remaining computations
    if(tid==0) {
        int remaining = data_pts%num_threads;
        for(int c1 = 0;c1<remaining;c1++) {
            double x1 = ((double)rand()) / ((double)RAND_MAX);
            double y1 = ((double)rand()) / ((double)RAND_MAX);
            double r1 = sqrt(x1*x1+y1*y1);
            if(r1<1) *circ_count += 1;
        }
    }
    
    //Now this value will be available in the join of pthread as status
    pthread_exit((void *)circ_count);
}
int main() {
    printf("Enter the number of threads(1-10): ");
    scanf("%d",&num_threads);
    printf("Enter the number of data points(10-1000000): ");
    scanf("%d",&data_pts);
    
    time_t t;
    //set the seed
    srand((unsigned) time(&t));
    
    //Create num_threads
    pthread_t threads[num_threads];
    //Return code
    int rc;
    //Return value from thread exit
    void *status;
    //total number of points inside the circle
    double tot_circle = 0;
    
    for(int i=0;i<num_threads;i++) {
        rc = pthread_create(&threads[i], NULL, func, (void *)(long)i);
        //Check if thread was created successfully
        if(rc) {
            printf("Error code: %d",rc);
            exit(-1);
        }
    }
    
    //Join the threads
    for(int i=0;i<num_threads;i++) {
        //status stores the final value calculated by the thread on thread_exit
        pthread_join(threads[i], &status);
        
        tot_circle += *(double *) status;
    }
    
    //Print the value of pi by monte carlo simulation
    printf("Value of Pi:%f",4*(tot_circle/data_pts));
    
    pthread_exit(NULL);
    return 0;
}

Compiling the program:

let the name of the file is monte.c

Then, to compile:

gcc -lm -pthread monte.c -o monte

To run:

./monte

Sample output:

Enter the number of threads(1-10): 10                                                                                         
Enter the number of data points(10-1000000): 1000000                                                                          
Value of Pi:3.142976 

Note: On some machines, C int is of 2 bytes. In that case use long int as it cannot store 1000000.