Calculating PI using Monte Carlo method gives imprecise answer

Question

I am trying to calculate PI using Monte Carlo method. My code gives the result 3.000 no matter how big MAXLEN is. After debugging it many times, I couldn't get what I'm doing wrong.

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

#define sqr2(x) ((x)*(x))
#define frand() ((double) rand() / (RAND_MAX))
#define MAXLEN 1000


int circumscribed(int radius){
    float xcoord = frand();
    float ycoord = frand(); 
    float coord = sqr2(xcoord) + sqr2(ycoord);

    if(coord <= radius)
        return 1;
    return -1;      
}

int main()
{
    int i;
    int circles = 0, rect = 0;;
    for(i = 0; i < MAXLEN; i++)
    {
        if(circumscribed(1) > 0)   // if(circumscribed(1)) shoul be enough but it doesn't work. Very odd in my opinion.
            circles++;
        rect++;  //this is equal to MAXLEN, I just used it for debugging
    }

    float PI = 4 * circles / rect;
    printf("PI is %2.4f: \n", PI);
    return 0;   
}

I'm not using C for years but try `float PI = 4.0 * circles / rect;`. — Jokester, Jul 15 '13 at 16:59
possible duplicate of [divide floats with integers in c](http://stackoverflow.com/questions/13331054/divide-floats-with-integers-in-c) — Oliver Charlesworth, Jul 15 '13 at 16:59
Can you say ["Integer Truncation"](http://en.wikipedia.org/wiki/Round-off_error) ;)? — paulsm4, Jul 15 '13 at 17:00
You might also want to consider seeding the random number generator before you use it. — John Sheridan, Jul 15 '13 at 17:04
Tthat solved it, thanks all. I'll accept the answer by @Bill after 7 minutes. — Abraham Guchi, Jul 15 '13 at 17:04

score 5 · Accepted Answer · answered Jul 15 '13 at 16:59

5

Since circles and rect are both int, the result of 4 * circles / rect will be an int. Use floating point numbers instead.

float PI = 4.0 * (float)circles / rect;

answered Jul 15 '13 at 16:59

Bill the Lizard

398,270
210
566
880

abelenky · Answer 2 · 2013-07-15T17:05:09.577

This expression is all integers:

4 * circles / rect;

Therefore, the result is an integer (3 in this case).

(as a similar example: 10 / 3 == 3, but 10.0 / 3.0 == 3.333333)

Try instead:

4.0 * circles / rect;

Just changing the (int)4 to a (double)4 by referring to it as 4.0 or even 4. should be enough.

Other Misc Observations

This line has an extra semi-colon:

int circles = 0, rect = 0;;

Your function circumscribed uses float. Your variable PI is also a float.
If you use double, you'll get greater precision.

score 2 · Answer 3 · answered Jul 15 '13 at 17:12

2

You are doing integer Math here. circles and rect are both int in your code, so the result of 4 * circles / rect is also an int. Use floating point numbers instead. Use double for better precision.

double PI = 4.0 * (double)circles / rect;

answered Jul 15 '13 at 17:12

Ankit Aggarwal

2,367
24
30

1

FYI, I've rolled back some of your previous edits as per these discussions: http://meta.stackexchange.com/q/135112/217057 http://meta.stackexchange.com/q/165700/217057 http://meta.stackexchange.com/q/112275/217057 – 000 Jul 15 '13 at 17:14
+1 for using `double` as 4.0 is a double and the printf() promotes the PI argument, if left as a float, to a double anyways. – chux - Reinstate Monica Jul 15 '13 at 18:49

score 1 · Answer 4 · answered Jul 15 '13 at 17:03

1

You are performing integer math here:

float PI = 4 * circles / rect;

Changing 4 to 4.0 is sufficient to fix the problem:

float PI = 4.0 * circles / rect;

Working version here

answered Jul 15 '13 at 17:03

Shafik Yaghmour

154,301
39
440
740

score 1 · Answer 5 · answered Jul 15 '13 at 19:39

float PI = 4 * circles / rect;

This performs integer math to the right of the = thus limiting your result to 1 significant digit.

Instead:

double PI = 4.0 * circles / rect;  // Best

Details

float PI = 4.0 * (float)circles / rect;  // OK

Recommend avoid using float. Use double instead to avoid limiting this Monte Carlo to about 7 digits. As circles becomes large, it may not convert exactly into the same float value. Instead it converts into a rounded float. This happens at about circles > 8,000,000 on many machines. By rounding, you are unnecessarily limiting the attainable precision. Using a double, this rounding does not occur until circles is about 9e15.

float PI = 4.0 * (double) circles / rect; // Better

The explicit cast of (double) circles may be useful to the reader, but code will perform the same way with or without it. 4.0 is a double and will cause a double promotion to circles before the multiplicity occurs, even without the cast.

double PI = 4.0 * circles / rect; // Better

As the 4.0 * ... result is a double, best precision is retained by saving to a double. Using float PI causes the multiplication/division, which was done in double precision to reduce its precision on saving to a float.

float or double PI = ...
printf("PI is %2.4f: \n", PI);

Note: Here, a PI converted to a double is passed to printf() regardless if PI was declared float or double. So might as access the precision in PI by declaring it double.

Note: int, float and double range and precision are machine dependent. The above reflects a common implementation.

score 0 · Answer 6 · answered Aug 02 '17 at 07:46

CHECK OUT THIS CODE:

#include <stdlib.h>
#include <stdio.h>
#include <math.h>
#include <string.h>
#define SEED 35791246

int main(int argc, char** argv)
{
   int niter=0;
   double x,y;
   int i,count=0; /* # of points in the 1st quadrant of unit circle */
   double z;
   double pi;

   printf("Enter the number of iterations used to estimate pi: ");
   scanf("%d",&niter);

   /* initialize random numbers */
   srand(SEED);
   count=0;
   for ( i=0; i<niter; i++) {
      x = (double)rand()/RAND_MAX;
      y = (double)rand()/RAND_MAX;
      z = x*x+y*y;
      if (z<=1) count++;
      }
   pi=(double)count/niter*4;
   printf("# of trials= %d , estimate of pi is %g \n",niter,pi);

return 0;
}

Calculating PI using Monte Carlo method gives imprecise answer

6 Answers6

Linked