Pthreads to compute Pi

[ndpndp]

I am planning to use parallelization for the code shown below to compute the value of a pi. Can anyone help me in guiding me to turn this code into multiple threads?

#include <pthread.h>
#include <semaphore.h>
#include <unistd.h>
#include <iostream>
#include <iomanip>
#include <time.h>
#include <cmath>
#include <stdlib.h> 
#define SHARED 1
#define MAX_THREADS     200
#define totalpoints 100000000
#endif

using namespace std; 
int main() 
{ 
   //Variables  
   double pi, 
   pointsIn = 0, 
   xValue, 
   yValue, 
   distance; 
   srand(time(NULL)); 
   for(int i=0; i<totalpoints; i++)
   { 
      xValue = (double) rand()/RAND_MAX; //Random point: x 
      yValue = (double) rand()/RAND_MAX; //Random point: y
      distance = (xValue * xValue) + (yValue * yValue); 
      distance = sqrt(distance); 
      if(distance <= 1) 
         pointsIn++; 
   } 
   pi = 4 * pointsIn/totalpoints; 
   cout<<"The vaule of pi according to monte carlo method is: " << pi << "."; 
   system("PAUSE");    
   return 0; 
} 

[kbw]

That's using the Monte Carlo method. Just do the same thing in any number of threads. When they're all done, sum pointsIn and totalpoints for each thread and do the final calculation pi = 4 * pointsIn/totalpoints at the end.

Monte Carlo converges very slowly.

[ndpndp]

Thanks for the reply, however i might not be clear in previous question. I want to divide the square into 4 parts and integrate at the end.

[kbw]

You only need to do the calculation on one quadrant and multiply by 4 at the end, that is the algorithm. If you're going to use threads, you just want to run more iterations in other threads.
http://math.fullerton.edu/mathews/n2003/montecarlopimod.html

Are you able to create a thread to do anything at all or are you stuck on threads.

I just noticed you're taking the square root to get the hypotenuse. The whole point of using a unit circle (circle with radius 1) is so you can avoid that expensive square root operation.

[ndpndp]

I am stuck on the threads and not able to move on, i want to use pthreads to run the program little faster.

[hanst99]

If you want to make this faster, before thinking about multithreading first think about how you're gonna parallelize your algorithm. There's no point in multithreading for solving a single problem if you can't separate the program into parts that can be calculated independently.

[ndpndp]

Now i am able to understand the problem a little bit. I want to calculate the values in 4 different threads and add them together to get a value of Pi. But i don't have any leads on how to do it.

[kev82]

You've already mentioned you want to use pthreads, have you considering googling for something like

"pthreads tutorial" or "pthreads montecarlo"

Took me roughly 30s to find:

http://myhomepage.ferris.edu/~nystroj/cpsc390/monte.cpp

[kbw]

First of all, I think you need to understand that Monte Carlo PI algorithm. I get the impression you don't quite get it.

Once you do, you'll see that there is no dependency between threads as far as the algorithm goes. The only thing each thread needs to calculate is the number of points in the the quadrant. You then find the ratio and multiply by 4 at the end in the main thread.

Here's an empty pthread example without your algorithm. It should be trivial to complete the task now.

#include <pthread.h>
#include <iostream>
#include <stdlib.h>

const size_t N = 4;

struct Info
{
	Info() : counter(0) {}

	unsigned counter;
};

void* func(void* param)
{
	if (Info* info = reinterpret_cast<Info*>(param))
	{
		// do some work
		info->counter = rand();
	}

	return param;
}

int main()
{
	srand(2012);

	pthread_t handles[N];
	Info info[N];

	// start the threads
	for (size_t i = 0; i != 4; ++i)
		pthread_create(&handles[i], NULL, func, &info[i]);

	// wait for the threads to complete
	for (size_t i = 0; i != 4; ++i)
		pthread_join(handles[i], NULL);

	// process the reply
	for (size_t i = 0; i != 4; ++i)
		std::cout << "value[" << i << "] = " << info[i].counter << std::endl;

	return 0;
}

[Cubbi]

if (Info* info = reinterpret_cast<Info*>(param))

A bit of nitpicking since the result is the same in this simple case, bit this is the wrong cast, static_cast<Info*> is the correct one (or C-style explicit cast, which resolves to static in this case)

[ndpndp]

Thanks for all the replies, however i am getting an error in the program given in http://myhomepage.ferris.edu/~nystroj/cpsc390/monte.cpp

it says something like:

cast from ‘void*’ to ‘int’ loses precision -fpermissive

I use g++ in linux, i googled out and changed it to intptr_t, but still the same error

[kbw]

You've gone completely overboard with that solution.

I can see where you use that PI constant in your program, but as your program is calculating PI, you probably shouldn't have one unless it's a requirement to print the actual precision.

Each of the threads can run the Monte Carlo algorithm concurrently, there is no interaction between the threads, so there's no reason to sync threads at all. IMHO, all those mutexes are surplus. You already stated what the algorithm is in your first post.

   for(int i=0; i<totalpoints; i++)
   { 
      double xValue = (double) rand()/RAND_MAX; //Random point: x 
      double yValue = (double) rand()/RAND_MAX; //Random point: y
      double distance = (xValue * xValue) + (yValue * yValue); 
      distance = sqrt(distance); // we can remove this expensive redundant line
      if (distance <= 1.0) 
         ++pointsIn;  // use pre-increment when you can
   }

I've added my tweaks in bold.

At the end, you add the pointsIn from each thread to get pointsInTotoal, then the final calculation is:

  pi = 4.0 * pointsInTotal/(NTHREADS * totalpoints);

Cubbi: Point taken.

[ndpndp]

Dear kbw,

I really appreciate your advice. As i am new to programming field, i am willing to implement the parallelization in my research for making my simulations faster. All i need at this moment is to run all 4 threads concurrently for parallelizing the problem and learn how to do it. I have eclipse running on ubuntu. May be my learning skills in programming are not so good at the moment.

[kbw]

I understand, we all have to start somewhere.

The problem really isn't as complicated as you've made it. You already have a template thread program where each thread can populate it's own element of an array of records. And you already started off with the algorithm that you need to use.

In most instances of multi-threaded programming, there are timing and data sharing issues that need to be maanged. Your timing issue is solved by the main thread starting all threads and waiting for them to complete before attempting to use the results. And you have no data sharing issues because the calculations can run completely independently of each other.

I think you're getting stuck because of the algorithm. Can you explain in your own words the idea of the algorithm? We can work on that, and once you understand it, it should all come together I think.

[ndpndp]

 If i use just this particular thing: double distance = (xValue * xValue) + (yValue * yValue);

We are not considering the hypotenuse as i am doing in the following case  
      distance = sqrt(distance);

Here in this algorithm if the distance is less than 1, it should be with in the circle else it is inside the square, finally we can divide the points inside the circle (pointsInTotal) and the total points (square). And able to estimate the value of Pi.

Please do correct me if i am wrong

My main problem is in implementation, in particular about the syntax

[kbw]

The idea is to use a game of chance to esitmate a value for PI.

You start by drawing a square of length 1 with a quarter of a circle of radius 1 inside the square.

If you were to throw a dart randomly at this shape, what are the chances of it landing inside the quarter circle?

The answer is (area of quarter circle) / (area of square). As the area of a square of length 1 is 1, the chance becomes (area of quarter circle) which is (PI * radius squared / 4). And as the radius of 1 squared is still 1, the answer is PI/4.

So the chance of your dart landing in the quarter circle is PI/4. Don't read further until you understand this point.

If we do this once, we'll get a value of 0 or 4 for PI. But as we do it more and more, we get better approximations for PI.

You could ask 3 friends to help you do this test be throwing darts and checking how many were in the quarter circle for a fixed number of throws. At the end, you'd find the ratio of how many were in, then multiply that by 4 and get a value for PI.

Back to the computer program. We throw the dart by getting random (x,y) coordinates. We know it's in the circle if the distance from the origin (0,0) is less than or equal to 1. And we calculate the distance using Pythagoras' theorem, c² = a² + b². But as 1² is still 1, we can just work c² and not attempt the square root.

Each thread does the same test and you just check the results at the end.

I hope you get from this:
1. how the algorithm works
2. see that the more throws that happen, the more accurate the result is
3. the throws can all happen independently/concurrently, you just need the ratio of in/total for each thrower
4. you do a simple calculation at the end, 4*sum of ratios to get PI
5. if you don't understand what to do, you can't instruct a person or machine to do it either

[ndpndp]

Thats a pretty good explanation, i finally got it and trying to change the code and it looks like:

#include <semaphore.h>
#include <unistd.h>
#include <iomanip>
#include <cmath>
#include <pthread.h>
#include <iostream>
#include <stdlib.h>

const size_t N = 4;

struct Info
{
	Info() : pi(0) {}

	unsigned pi;
};

void* func(void* param)
{
	double pi, pointsIn = 0,
			   xValue,
			   yValue,
			   totalpoints=1000000,
			   distance;
			   srand(time(NULL));

			   if (Info* info = static_cast<Info*>(param))
	{


		   //Monte-carlo method to calculate the value of Pi
		   for(int info=0; info<totalpoints; info++)
		   {
		      xValue = (double) rand()/RAND_MAX; //Random point: x
		      yValue = (double) rand()/RAND_MAX; //Random point: y
		      distance = (xValue * xValue) + (yValue * yValue);
		      if(distance <= 1.0)
		         ++pointsIn;
		   }
		   pi = 4.0 * pointsInTotal/(4 * totalpoints);
	}

	return param;
}

int main()
{
	srand(time(NULL));

	pthread_t handles[N];
	Info info[N];

	// start the threads
	for (size_t i = 0; i != 4; ++i)
		pthread_create(&handles[i], NULL, func, &info[i]);

	// wait for the threads to complete
	for (size_t i = 0; i != 4; ++i)
		pthread_join(handles[i], NULL);

	// process the reply
	for (size_t i = 0; i != 4; ++i)
		std::cout << "value[" << i << "] = " << info[i].pi << std::endl;

	return 0;
}

[kbw]

The algorithm as I described it is for each thrower to throw they're darts independently, then the totals a brought together and PI calculated from the results of all throwers. So in total the calculation is based on (N*totalpoints) throws.

Your code makes each thrower attempt to calculate PI from his own results. So in total your calculation is based on (totalpoints) throws, with N independent results. It's not quite the same thing.

Once all the threads complete, you should get the totals from all threads and calculate PI from that.

[ndpndp]

Actually i am getting the result as:

value[0] = 0
value[1] = 0
value[2] = 0
value[3] = 0

[kbw]

That's because you're not returning the calculation result, but that isn't really the problem. Has it not occurred to you that you should only have one value for PI, and not 4?