Parallel program to calculate PI

[Prime Mover]

Hello all,

I have got the pseudo-code below that I would like to convert to c
language. The algorithm calculates Pi value. I am somewhat familiar
with C language, but I am just starting to learn parallel programming.
In this specific example, the idea seems to be simple: one must
subdivide the main loop into pieces that can be executed by
independent "tasks" (computers??).

Then, each "worker task" executes a part of the loop a certain number
of times, independently of the other worker tasks. One specific task
plays the role of "master task", which will collect and sum the
results of the worker tasks:

% descriptive algorithm:
1. Inscribe a circle inside a square
2. Generate random points inside the square
3. Determine the number of points that fell inside the circle
4. Let r be the number of points inside the circle divided by the
total number of points
5. Pi is approximately equal to 4*r
6. The more points are generated, the more is the precision in P value

% pseudo-code (parallel):
1. npoints = 10000
2. circle_count = 0
3. p = number of tasks
4. num = npoints/p
5. find out if I am MASTER or WORKER
6. do j = 1,num
7. generate 2 random numbers between 0 and 1
8. xcoordinate = random1
9. ycoordinate = random2
10. if (xcoordinate, ycoordinate) inside circle then
circle_count = circle_count + 1
11. end do
12. if I am MASTER
13. receive from WORKERS their circle_counts
14. compute PI (use MASTER and WORKER calculations)
15. else if I am WORKER
16. send to MASTER circle_count
17. endif

Any help would be much appreciated.

I thank you all in advance.

[Prime Mover]

Let me just be a bit more specific:

My (understading) problem starts in the line 5 of the pseudo-code:
5. find out if I am MASTER or WORKER

How would I specificy a "worker"? Would that be another computer?
If yes, how can I access this remote computer in the calculations, in
C language?
If yes, it means that I have to have a LAN or something to perform
tests?

I have found that there are some libraries such as OMP or MPI that
could be
used, but I'd like to know if there is a more "raw" way of doing this
first.

Again, thank you all.

[Spiros Bousbouras]

On 13 May, 15:20, Prime Mover <[e-mail]:

Hello all,

I have got the pseudo-code below that I would like to convert to c
language. The algorithm calculates Pi value. I am somewhat familiar
with C language, but I am just starting to learn parallel programming.
In this specific example, the idea seems to be simple: one must
subdivide the main loop into pieces that can be executed by
independent "tasks" (computers??).

Then, each "worker task" executes a part of the loop a certain number
of times, independently of the other worker tasks. One specific task
plays the role of "master task", which will collect and sum the
results of the worker tasks:

% descriptive algorithm:
1. Inscribe a circle inside a square
2. Generate random points inside the square
3. Determine the number of points that fell inside the circle
4. Let r be the number of points inside the circle divided by the
total number of points
5. Pi is approximately equal to 4*r
6. The more points are generated, the more is the precision in P value

% pseudo-code (parallel):
1. npoints = 10000
2. circle_count = 0
3. p = number of tasks
4. num = npoints/p
5. find out if I am MASTER or WORKER
6. do j = 1,num
7. generate 2 random numbers between 0 and 1
8. xcoordinate = random1
9. ycoordinate = random2
10. if (xcoordinate, ycoordinate) inside circle then
circle_count = circle_count + 1
11. end do
12. if I am MASTER
13. receive from WORKERS their circle_counts
14. compute PI (use MASTER and WORKER calculations)
15. else if I am WORKER
16. send to MASTER circle_count
17. endif

On 13 May, 15:28, Prime Mover <[e-mail]:

Let me just be a bit more specific:

My (understading) problem starts in the line 5 of the pseudo-code:
5. find out if I am MASTER or WORKER

How would I specificy a "worker"? Would that be another computer?
If yes, how can I access this remote computer in the calculations, in
C language?
If yes, it means that I have to have a LAN or something to perform
tests?

I have found that there are some libraries such as OMP or MPI that
could be
used, but I'd like to know if there is a more "raw" way of doing this
first.

Standard C has no built-in support for parallel processing so an
answer
to the question "find out if I am MASTER or WORKER" falls outside
standard C which is what's topical here. I have no idea what kind of
hardware set-up and extensions to C can be used to tackle numerically
intensive parallel algorithms. Perhaps others here have the relevant
experience. Searching Google groups for *parallel* I found
comp.parallel.mpi and comp.parallel.pvm where I'm guessing you might
get more useful advice than here.

Out of curiosity for what value of npoints are you aiming for ? In
your example it's only 100000 and you can get that on a modern desktop
in a few seconds.

[Spiros Bousbouras]

On 13 May, 15:28, Prime Mover <[e-mail]:

>

Let me just be a bit more specific:

My (understading) problem starts in the line 5 of the pseudo-code:
5. find out if I am MASTER or WORKER

How would I specificy a "worker"? Would that be another computer?
If yes, how can I access this remote computer in the calculations, in
C language?
If yes, it means that I have to have a LAN or something to perform
tests?

I have found that there are some libraries such as OMP or MPI that
could be
used, but I'd like to know if there is a more "raw" way of doing this
first.

Standard C has no built-in support for parallel processing so an
answer
to the question "find out if I am MASTER or WORKER" falls outside
standard C which is what's topical here. I have no idea what kind of
hardware set-up and extensions to C can be used to tackle numerically
intensive parallel algorithms. Perhaps others here have the relevant
experience. Searching Google groups for *parallel* I found
comp.parallel.mpi and comp.parallel.pvm where I'm guessing you might
get more useful advice than here.

There's also comp.parallel

[Pietro Cerutti]

Prime Mover wrote:

% descriptive algorithm:
1. Inscribe a circle inside a square
2. Generate random points inside the square
3. Determine the number of points that fell inside the circle

How can you determine it without knowing PI a priori?

Pietro Cerutti

[Bart]

On May 13, 3:55*pm, Spiros Bousbouras <[e-mail]:

On 13 May, 15:20, Prime Mover <[e-mail]:

I have got the pseudo-code below that I would like to convert to c
language. The algorithm calculates Pi value. I am somewhat familiar

Out of curiosity for what value of npoints are you aiming for ? In
your example it's only 100000 and you can get that on a modern desktop
in a few seconds.- Hide quoted text -

I got 3.1415 using a billion points. Looks like it will converge very
slowly.

Also the granularity of the x,y points may affect the maximum accuracy
(because it introduces errors near the circular edge). Tried a sphere
too but not any better.

--
Bartc

[Richard Heathfield]

Pietro Cerutti said:

Prime Mover wrote:

>% descriptive algorithm:
1. Inscribe a circle inside a square
2. Generate random points inside the square
3. Determine the number of points that fell inside the circle

How can you determine it without knowing PI a priori?

There's nothing in the above description that requires a knowledge of pi.

For example, you can draw a circle of radius r easily enough just by
applying the equation x^2 + y^2 = r^2 (no, that is not C code!) for all y
in the range -r to +r, taking the positive result and its negative
analogue for your x points. (Actually, it's even easier than that, because
there are other symmetries to exploit.)

Determining whether you're in the circle can be done either by tedious
inspection (a process not unlike floodfill) or (much faster) by
determining whether x^2 + y^2 <= r^2.

--
Richard Heathfield <http://www.cpax.org.uk>
Email: -http://www. +rjh@
Google users: <http://www.cpax.org.uk/prg/writings/googly.php>
"Usenet is a strange place" - dmr 29 July 1999

[Spiros Bousbouras]

On 13 May, 16:02, Pietro Cerutti <gahr@localhostwrote:

Prime Mover wrote:

% descriptive algorithm:
1. Inscribe a circle inside a square
2. Generate random points inside the square
3. Determine the number of points that fell inside the circle

How can you determine it without knowing PI a priori?

Imagine a circle with radius is 1 and its center has coordinates
(0 , 0).
A random point with coordinates (x,y) will be inside the square
if -1 <= x <= 1 and -1 <= y <= 1
if ( x*x + y*y < 1) /* inside the circle */

[Spiros Bousbouras]

On 13 May, 16:15, Bart <[e-mail]:

On May 13, 3:55 pm, Spiros Bousbouras <[e-mail]:

On 13 May, 15:20, Prime Mover <[e-mail]:

I have got the pseudo-code below that I would like to convert to c
language. The algorithm calculates Pi value. I am somewhat familiar

Out of curiosity for what value of npoints are you aiming for ? In
your example it's only 100000 and you can get that on a modern desktop
in a few seconds.- Hide quoted text -

I got 3.1415 using a billion points. Looks like it will converge very
slowly.

Also the granularity of the x,y points may affect the maximum accuracy
(because it introduces errors near the circular edge). Tried a sphere
too but not any better.

And of course you must know that your random numbers
will be uniformly distributed within the square.

[CBFalconer]

Pietro Cerutti wrote:

Prime Mover wrote:

>% descriptive algorithm:
1. Inscribe a circle inside a square
2. Generate random points inside the square
3. Determine the number of points that fell inside the circle

How can you determine it without knowing PI a priori?

Write the expression for the area of a circle. Do the same for the
area of the superscribed square. Take the ratio of areas. It
contains PI.

--
[mail]: Chuck F (cbfalconer at maineline dot net)
[page]: <http://cbfalconer.home.att.net>
Try the download section.

** Posted from http://www.teranews.com **

[Richard Heathfield]

Spiros Bousbouras said:

<snip>

And of course you must know that your random numbers
will be uniformly distributed within the square.

That's easy. Use the following random point generator:

#include <assert.h>

void rndpt(unsigned long int *x,
unsigned long int *y,
unsigned long int max) /* max = side of square */
{
static unsigned long int n = 0;
assert(x != NULL && y != NULL);
*x = n % max;
*y = n++ / max;
n %= max;
return;
}

If you call this i * max times, where max is a constant and i is an
integer, the distribution of the random points will be uniform.

--
Richard Heathfield <http://www.cpax.org.uk>
Email: -http://www. +rjh@
Google users: <http://www.cpax.org.uk/prg/writings/googly.php>
"Usenet is a strange place" - dmr 29 July 1999

[Spiros Bousbouras]

On 13 May, 17:00, Richard Heathfield <[e-mail]:

Spiros Bousbouras said:

<snip>

And of course you must know that your random numbers
will be uniformly distributed within the square.

That's easy. Use the following random point generator:

#include <assert.h>

void rndpt(unsigned long int *x,
unsigned long int *y,
unsigned long int max) /* max = side of square */
{
static unsigned long int n = 0;
assert(x != NULL && y != NULL);
*x = n % max;
*y = n++ / max;
n %= max;
return;

}

If you call this i * max times, where max is a constant and i is an
integer, the distribution of the random points will be uniform.

Since *y will always get the value 0 I don't think
so.

[Richard Heathfield]

Spiros Bousbouras said:

On 13 May, 17:00, Richard Heathfield <[e-mail]:

>Spiros Bousbouras said:

<snip>

And of course you must know that your random numbers
will be uniformly distributed within the square.

That's easy. Use the following random point generator:

#include <assert.h>

void rndpt(unsigned long int *x,
unsigned long int *y,
unsigned long int max) /* max = side of square */
{
static unsigned long int n = 0;
assert(x != NULL && y != NULL);
*x = n % max;
*y = n++ / max;
n %= max;
return;

}

If you call this i * max times, where max is a constant and i is an
integer, the distribution of the random points will be uniform.

Since *y will always get the value 0 I don't think
so.

Since it won't, I do.

--
Richard Heathfield <http://www.cpax.org.uk>
Email: -http://www. +rjh@
Google users: <http://www.cpax.org.uk/prg/writings/googly.php>
"Usenet is a strange place" - dmr 29 July 1999

[Eric Sosman]

Richard Heathfield wrote:

Spiros Bousbouras said:

>On 13 May, 17:00, Richard Heathfield <[e-mail]:

>>Spiros Bousbouras said:

<snip>

And of course you must know that your random numbers
will be uniformly distributed within the square.
That's easy. Use the following random point generator:

#include <assert.h>

void rndpt(unsigned long int *x,
unsigned long int *y,
unsigned long int max) /* max = side of square */
{
static unsigned long int n = 0;
assert(x != NULL && y != NULL);
*x = n % max;
*y = n++ / max;
n %= max;
return;

}

If you call this i * max times, where max is a constant and i is an
integer, the distribution of the random points will be uniform.

Since *y will always get the value 0 I don't think
so.

Since it won't, I do.

Try running it for a while and making a histogram of
the y outputs.

(The bug is in the line preceding the return.)

--
Er*********@sun.com

[Bart]

On May 13, 5:00*pm, Richard Heathfield <[e-mail]:

Spiros Bousbouras said:

<snip>

And of course you must know that your random numbers
will be uniformly distributed within the square.

That's easy. Use the following random point generator:

#include <assert.h>

void rndpt(unsigned long int *x,
* * * * * *unsigned long int *y,
* * * * * *unsigned long int max) /* max = side of square */
{
* static unsigned long int n = 0;
* assert(x != NULL && y != NULL);
* *x = n % max;
* *y = n++ / max;
* n %= max;
* return;

}

If you call this i * max times, where max is a constant and i is an
integer, the distribution of the random points will be uniform.

You mean max*max?

This looks to be simply filling in a square sequentially. You are then
effectively calculating the area of a circle in a square by counting
the all dots. That's not really in the spirit of the original method.
(I think it also converges more slowly compared with the same number
of points picked randomly.)

--
Bartc

[dj3vande]

In article <5s******************************@bt.com>,
Richard Heathfield <rj*@see.sig.invalidwrote:

>Spiros Bousbouras said:

>On 13 May, 17:00, Richard Heathfield <[e-mail]:

>>Spiros Bousbouras said:

<snip>

And of course you must know that your random numbers
will be uniformly distributed within the square.

That's easy. Use the following random point generator:

#include <assert.h>

void rndpt(unsigned long int *x,
unsigned long int *y,
unsigned long int max) /* max = side of square */
{
static unsigned long int n = 0;
assert(x != NULL && y != NULL);
*x = n % max;
*y = n++ / max;
n %= max;
return;

}

If you call this i * max times, where max is a constant and i is an
integer, the distribution of the random points will be uniform.

Since *y will always get the value 0 I don't think
so.

Since it won't, I do.

It looks to me like it will, so I don't.

When n is max-1, the line that assigns *y sets *y to 0, and then bumps
n to max. The next line folds max back to 0, so on the next invocation
you'll go back to 0 for *x and *y will still be 0.
I think you meant to write "n %= (max*max)" for the line just before
the return.
dave

--
Dave Vandervies dj3vande at eskimo dot com

I'd like to believe that somewhere there must be a BMW not driven by a
fsckwit. --Garrett Wollman and Richard P. Grant in

Give me a BMW and I'll fulfill that wish. the scary devil monastery

[Richard Heathfield]

Eric Sosman said:

Richard Heathfield wrote:

>Spiros Bousbouras said:

>>On 13 May, 17:00, Richard Heathfield <[e-mail]:

<snip>

>>> n %= max;
return;

}

If you call this i * max times, where max is a constant and i is an
integer, the distribution of the random points will be uniform.
Since *y will always get the value 0 I don't think
so.

Since it won't, I do.

Try running it for a while and making a histogram of
the y outputs.

(The bug is in the line preceding the return.)

My apologies to Spiros - he's right and I'm wrong. It should have been:

n %= max * max;

Furthermore, it has to be called i * max * max times, not i * max times,
for the distribution to be uniform.

--
Richard Heathfield <http://www.cpax.org.uk>
Email: -http://www. +rjh@
Google users: <http://www.cpax.org.uk/prg/writings/googly.php>
"Usenet is a strange place" - dmr 29 July 1999

[Richard Heathfield]

Bart said:

On May 13, 5:00 pm, Richard Heathfield <[e-mail]:

<snip>

>n %= max

* max;

>;
return;

}

If you call this i * max times, where max is a constant and i is an
integer, the distribution of the random points will be uniform.

You mean max*max?

Yes.

This looks to be simply filling in a square sequentially.

Yes.

You are then
effectively calculating the area of a circle in a square by counting
the all dots.

Yes.

That's not really in the spirit of the original method.

No - but at least it would have been funny, if only I hadn't screwed it up.

I once saw a wonderful pi calculation method on the Web (to which I can't
find a URL right now) - it was a Monte Carlo method, which used millions
of digits of pi as a PRNG, and concluded... that pi is about 3. Fabulous!

--
Richard Heathfield <http://www.cpax.org.uk>
Email: -http://www. +rjh@
Google users: <http://www.cpax.org.uk/prg/writings/googly.php>
"Usenet is a strange place" - dmr 29 July 1999

[user923005]

On May 13, 10:06*am, Richard Heathfield <[e-mail]:
[snip]

I once saw a wonderful pi calculation method on the Web (to which I can't
find a URL right now) - it was a Monte Carlo method, which used millions
of digits of pi as a PRNG, and concluded... that pi is about 3. Fabulous!

You can also calculate pi by dropping needles on a grid and seeing if
they touch the lines.
http://www.angelfire.com/wa/hurben/buff.html

I have calculated pi to several thousand digits by numerical
integration of arctangent identities.

However, I used C++ so it is not even tangentially topical.

So let's revive the Dik Winter pi programs:

The minimalistic C source code to calculate pi to 32372 digits reads:

/* Calculation of pi to 32372 decimal digits */
/* Size of program: 152 characters */
/* After Dik T. Winter, CWI Amsterdam */
unsigned a=1e4,b,c=113316,d,e,f[113316],g,h,i;
main(){for(;b=c,c-=14;i=printf("%04d",e+d/a),e=d%a)
while(g=--b*2)d=h*b+a*(i?f[b]:a/5),h=d/--g,f[b]=d-g*h;}

An even shorter version creates the first 16276 digits of pi:

/* Calculation of pi to 16276 decimal digits */
/* Size of program: 143 characters */
/* After Dik T. Winter, CWI Amsterdam */
int a=1e4,b,c=56980,d,e,f[56980],g,h,i;
main(){for(;b=c,c-=14;i=printf("%04d",e+d/a),e=d%a)
while(g=--b*2)d=h*b+a*(i?f[b]:a/5),h=d/--g,f[b]=d%g;}

[Prime Mover]

Very interesting discussions!

Have anyone used header files such as omp.h or mpi.h that seem
to call functions for parallel computing in C?

I promise this would be the last question about parallel programming
here :-)

[Prime Mover]

On 13 maio, 16:30, Prime Mover <[e-mail]:

Very interesting discussions!

Have anyone used header files such as omp.h or mpi.h that seem
to call functions for parallel computing in C?

I promise this would be the last question about parallel programming
here :-)

I forgot to ask if anyone could share these files with me, because I
couldn't find them anywhere for a direct download.

Thank you again.

[Keith Thompson]

Prime Mover <ep*****@hotmail.comwrites:

Very interesting discussions!

Have anyone used header files such as omp.h or mpi.h that seem
to call functions for parallel computing in C?

I promise this would be the last question about parallel programming
here :-)

Assuming that omp.h is for OpenMP, see <http://openmp.org/>.

For MPI, see comp.parallel.mpi.

--
Keith Thompson (The_Other_Keith) ks***@mib.org <http://www.ghoti.net/~kst>
Nokia
"We must do something. This is something. Therefore, we must do this."
-- Antony Jay and Jonathan Lynn, "Yes Minister"

[Eric Sosman]

Prime Mover wrote:

On 13 maio, 16:30, Prime Mover <[e-mail]:

>Very interesting discussions!

Have anyone used header files such as omp.h or mpi.h that seem
to call functions for parallel computing in C?

I promise this would be the last question about parallel programming
here :-)

I forgot to ask if anyone could share these files with me, because I
couldn't find them anywhere for a direct download.

This is Question 10.11 in the comp.lang.c Frequently
Asked Questions (FAQ) list, <http://www.c-faq.com/>, and
you're not going to like the answer.

--
Er*********@sun.com

[CBFalconer]

Richard Heathfield wrote:

>

.... snip ...

>
No - but at least it would have been funny, if only I hadn't screwed it up.

I once saw a wonderful pi calculation method on the Web (to which I can't
find a URL right now) - it was a Monte Carlo method, which used millions
of digits of pi as a PRNG, and concluded... that pi is about 3. Fabulous!

However, rounded to the nearest integer, it was absolutely
correct. :-)

--
[mail]: Chuck F (cbfalconer at maineline dot net)
[page]: <http://cbfalconer.home.att.net>
Try the download section.
** Posted from http://www.teranews.com **

[Richard Heathfield]

user923005 said:

On May 13, 10:06 am, Richard Heathfield <[e-mail]:
[snip]

>I once saw a wonderful pi calculation method on the Web (to which I
can't find a URL right now) - it was a Monte Carlo method, which used
millions of digits of pi as a PRNG, and concluded... that pi is about 3.
Fabulous!

You can also calculate pi by dropping needles on a grid and seeing if
they touch the lines.

Yes, I know. I once demonstrated this to my kids, using the tiles on our
kitchen floor as the grid. I just asked one of them what value we came up
with, and he'd completely forgotten not only the value but the entire
demonstration. So it's time to do it all over again! :-)

<snip>

So let's revive the Dik Winter pi programs:

Please don't, or I'll turn green.

--
Richard Heathfield <http://www.cpax.org.uk>
Email: -http://www. +rjh@
Google users: <http://www.cpax.org.uk/prg/writings/googly.php>
"Usenet is a strange place" - dmr 29 July 1999

[user923005]

On May 13, 2:57*pm, Richard Heathfield <[e-mail]:

user923005 said:

On May 13, 10:06 am, Richard Heathfield <[e-mail]:
[snip]

I once saw a wonderful pi calculation method on the Web (to which I
can't find a URL right now) - it was a Monte Carlo method, which used
millions of digits of pi as a PRNG, and concluded... that pi is about 3..
Fabulous!

You can also calculate pi by dropping needles on a grid and seeing if
they touch the lines.

Yes, I know. I once demonstrated this to my kids, using the tiles on our
kitchen floor as the grid. I just asked one of them what value we came up
with, and he'd completely forgotten not only the value but the entire
demonstration. So it's time to do it all over again! :-)

<snip>

So let's revive the Dik Winter pi programs:

Please don't, or I'll turn green.

Here is a version with less UB:

#include <stdio.h>
static unsigned a = 10000,
b,
c = 113316,
d,
e,
f[113316],
g,
h,
i;
int main(void)
{
for (; b = c, c -= 14; i = (unsigned) printf("%04u", e + d / a), e
= d % a)
while ((g = --b * 2))
d = h * b + a * (i ? f[b] : a / 5), h = d / --g, f[b] = d
- g * h;
return 0;
}

[dSpam]

On 13 Mai, 16:28, Prime Mover <[e-mail]:

My (understading) problem starts in the line 5 of the pseudo-code:
5. find out if I am MASTER or WORKER

Read the description of fork() (the absence of which in the C standard
documents lets you conclude that it's not part of the C standard
library) or ask in a group dedicated to your computing environment.