Count Pi in Streams

Question

I have two implementations of counting pi with the Monte Carlo method : with and without threads. A threadless implementation works very well, but the threading method has problems with accuracy and performance. Here is the code:

Without threads:

#include <stdio.h> #include <stdlib.h> #include <time.h> int main() { srand(time(NULL)); unsigned long N = 0, Nin = 0; float x,y; while(N < 2E+9) { x = rand()/((float)RAND_MAX + 1.0)*10.0 - 5.0; y = rand()/((float)RAND_MAX + 1.0)*10.0 - 5.0; if(x*x + y*y < 25.0) Nin += 1; N++; } long double pi = 4.0 * (long double)Nin / (long double)N; printf("\tPi1: %.20Lf\n\t%lu %lu\n", pi, Nin, N); return 0; } 

And with threads:

 #include <stdio.h> #include <stdlib.h> #include <time.h> #include <pthread.h> typedef struct { unsigned long Nin; unsigned long N; } nums; void pi_counter(nums* a) { float x,y; unsigned int N = 0, Nin = 0; while(N < 1E+9) { x = rand()/((float)RAND_MAX + 1.0)*10.0 - 5.0; y = rand()/((float)RAND_MAX + 1.0)*10.0 - 5.0; if(x*x + y*y < 25.0) Nin++; N++; } a -> Nin += Nin; a -> N += N; } int main() { pthread_t thread1, thread2, thread3; nums a; srand(time(NULL)); pthread_create( &thread1, NULL, pi_counter, &a ); pthread_create( &thread2, NULL, pi_counter, &a ); pthread_join( thread1, NULL ); pthread_join( thread2, NULL ); long double pi = 4.0 * (long double)a.Nin / (long double)aN; printf("\tPi2: %.20Lf\n\t%lu %lu\n", pi, a.Nin, aN); return 0; } 

Results:

 $ time ./pi2 Pi2: 3.14147154999999999995 1570735775 2000000000 real 1m1.927s user 1m23.624s sys 0m0.139s $ time ./pi Pi1: 3.14158868600000000006 1570794343 2000000000 real 0m49.956s user 0m49.887s sys 0m0.022s 

Where is my mistake?

Answer 1

rand not thread safe; simultaneous use in multiple threads will result in undefined behavior . You can either wrap it with a function that acquires and holds the mutex when calling rand , or you can use rand_r or (even better) write a decent PRNG to use in its place.

Answer 2

Among other answers, in the following code

 a -> Nin += Nin; a -> N += N; 

a is a common, but not protected mutex, which leads to improper addition. Although you may not have encountered this problem, but in the end.

Answer 3

Your rand() at the same time in the threads will give a sequence of the same numbers , because you get different results, but the algorithm is excellent (it is probabilistic , nothing is guaranteed). Why is the same sequence? Because an instance of a seed term is for each process, a thread is a process, but easy.