Estimating the value of Pi using the Monte Carlo Method

Question

I am writing a C program that will be able to accept an input value that dictates the number of iterations that will be used to estimate Pi.

For example, the number of points to be created as the number of iterations increases and the value of Pi also.

Here is the code I have so far:

#include <stdio.h>
#include <stdlib.h>
main()
{
    const double pp = (double)RAND_MAX * RAND_MAX;
    int innerPoint = 0, i, count;
    
    printf("Enter the number of points:");
    scanf("%d", &innerPoint);
    for (i = 0; i < count; ++i){
        float x = rand();
        float y = rand();
        
        if (x * x + y * y <= 1){
            ++innerPoint;
        }
        
        int ratio = 4 *(innerPoint/ i);
        
        printf("Pi value is:", ratio);
    }
}

Help fix my code as I'm facing program errors.

`rand` doesn't return a float, it returns an integer. See [this question](https://stackoverflow.com/questions/13408990/how-to-generate-random-float-number-in-c). — 1201ProgramAlarm, Jun 14 '21 at 20:53
Also, there is a [c-faq](http://c-faq.com/lib/rand48.html) entry. Also, consider fixed point? — Neil, Jun 14 '21 at 21:16
You may want to input the value of `count`, not `innerPoint` and compare x^2 + y^2 to `pp` (that *you* have introduced), not `1`. — Bob__, Jun 14 '21 at 21:17
`innerPoint/ i` is integer division, and `printf("Pi value is:", ratio);` does not output any value, whether or not it is correct. — Weather Vane, Jun 14 '21 at 22:04

score 1 · Answer 1 · answered Jun 14 '21 at 22:20

1

rand() returns an integer [0...RAND_MAX].

So something like:

   float x = rand()*scale;  // Scale is about 1.0/RAND_MAX

The quality of the Monte Carlo method is dependent on a good random number generator. rand() may not be that good, but let us assume it is a fair random number generator for this purpose.

The range of [0...RAND_MAX] is RAND_MAX+1 different values that should be distributed evenly from [0.0...1.0].

((float) rand())/RAND_MAX biases the end points 0.0 and 1.0 giving them twice the weight of others.

Consider instead [0.5, 1.5, 2.5, ... RAND_MAX + 0.5]/(RAND_MAX + 1).

RAND_MAX may exceed the precision of float so converting rand() or RAND_MAX, both int, to float can incurring rounding and further disturb the Monte Carlo method. Consider double.

#define RAND_MAX_P1 ((double)RAND_MAX + 1.0)

// float x = rand();
double x = ((double) rand() + 0.5)/RAND_MAX_P1;

x * x + y * y can also incur excessive rounding. C has hypot(x,y) for a better precision sqrt(x*x + y*y). Yet here, with small count, it likely makes no observable difference.

//  if (x * x + y * y <= 1)
if (hypot(x, y <= 1.0))

answered Jun 14 '21 at 22:20

chux - Reinstate Monica

143,097
13
135
256

+1 for correctly pointing out that `RAND_MAX` is inclusive and this requires `[0, 1)`. `double` is more precise; another way is to use only integers by scaling the problem up. – Neil Jun 14 '21 at 23:41
@neal True. An integer only soliton possible. – chux - Reinstate Monica Jun 15 '21 at 00:00
I implemented in `unsigned long int`-only; I think this gives 11-bit more accurate approximation close to magnitude 1. It's hard to say whether that is the bottleneck to achieving a more accurate `pi`. – Neil Jun 15 '21 at 08:48
An integer solution would be better if it is faster (more iterations in the same time should give a better estimate), but one has to be careful with overflow and possible non-uniform distribution if trying to do something like `rand() % MY_MAX_RAND` since `RAND_MAX` only comes with the guarantee that it is at least 32767. – nielsen Jun 15 '21 at 09:36
@nielsen If `a,b=rand();` then `x^2 + y^2 = (2a+1)^2 + (2b+1)^2 <= (2(RAND_MAX+1))^2` probably overflows, but you could cancel 8 and get the results, assuming `sizeof(long) >= 2*log RAND_MAX`, as a `long` exactly: `a * (a + 1) / 2 + b * (b + 1) / 2 <= (2ul + RAND_MAX) * RAND_MAX / 2` (`+3/4` dropped). – Neil Jun 15 '21 at 19:21
1

@Neil If `sizeof(long)>=2*sizeof(int)` then we should be ok using `unsigned long` for the calculations since `rand()` returns a non-negative signed `int`, but any solution assuming `sizeof(long)>sizeof(int)` would not be portable since the standard does not guarantee that. – nielsen Jun 16 '21 at 08:36

JANO · Answer 2 · 2021-06-14T22:09:17.870

I am sure it is not the best solution, but it should do the job and is similar to your code. Use a sample size of at least 10000 to get a value near PI.

As mentioned in the commenter: You should look at the data types of the return values functions give you.

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

int main()
{
    
    // Initialize random number generation
    srand(time(NULL));
    
    int samples = 10000;
    int points_inside =0;

    // Read integer - sample size (use at least 10000 to get near PI)
    printf("Enter the number of points:");
    scanf("%d", &samples);
    
    
    
    for (int i = 0; i < samples; ++i){
        
        // Get two numbers between 0 and 1
        float x = (float) rand() / (float)RAND_MAX;
        float y = (float) rand() / (float)RAND_MAX;
        
        
        // Check if point is inside
        if (x * x + y * y <= 1){
            points_inside++;
        }
                
        // Calculate current ratio
        float current_ratio = 4 * ((float) points_inside / (float) i);
        
        printf("Current value of pi value is: %f \n", current_ratio);
    }
}

Estimating the value of Pi using the Monte Carlo Method

2 Answers2