C: Approxmation of Pi with stratified sampling #include&lt;stdio.h&gt;
#include&lt;stdlib.h&gt;
#include&lt;time.h&gt;
#include&lt;math.h&gt;

int m, inbound, g;
double pi;
double uncertainty;
double pi_exactly=3.14159265358979323846;

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

    printf("Quantity of random points per square: \n");
    scanf("%d", &amp;m);

    printf("Gridsize: \n");
    scanf("%d", &amp;g);

    for(int i=0; i&lt;g; i++)
    {
        for(int j=0; j&lt;g; ++j)
        {

            if(pow((i+1) / g,2)+pow((j+1) / g,2) &lt; 1){
                inbound+=m;
            }else if(pow(i,2)+pow(j,2)&gt;1){
                break;
            }else{

                for(int k=1; k&lt;=m; k++)
                {
                    double x = ((double)rand() / RAND_MAX + i) * 1 / (double) g;
                    double y = ((double)rand() / RAND_MAX + j) * 1 / (double) g;

                    if(pow(x,2)+pow(y,2)&lt;=1){
                        inbound+=k;
                    }
                }
            }
        }
    }


    pi = (double) inbound / ((double) m *  pow(g,2))*4;
    uncertainty = pi - pi_exactly;

    printf("The approximation of pi is: %lf\n", pi);
    printf("Uncertainty: %e\n", uncertainty);
return 0;

}


I want to approximate \$\pi\$ with the Monte-Carlo-Simulation, so if the random points are inside the unit circle, then there are "inbound". Otherwise not. Pi can be randomize with the divison of the inbound points and the whole amount of random points. Unfortunately the code does not approximate \$\pi\$ very precisely. With \$m=100, g=1000\$ I get \$\pi = 3.991...\$ Is there something to improve? Thanks!
Bug #1

In your first test:


        if(pow((i+1) / g,2)+pow((j+1) / g,2) &lt; 1){



you didn't convert (i+1)/g to a double correctly.  You are using integer math, so (i+1)/g is going to be zero until i reaches g-1.  What you wanted to do was this:

        if(pow((double)(i+1) / g,2)+pow((double)(j+1) / g,2) &lt; 1){


Bug #2

In your second test:


        }else if(pow(i,2)+pow(j,2)&gt;1){



you forgot to divide by g.  So it should have been:

        }else if(pow((double)i/g,2)+pow((double)j/g,2)&gt;1){


Bug #3

Your last test was OK, but you added the wrong amount to inbound:


                if(pow(x,2)+pow(y,2)&lt;=1){
                    inbound+=k;
                }



Here, it should have been:

                if(pow(x,2)+pow(y,2)&lt;=1){
                    inbound++;
                }


Otherwise, later when you divide by the number of samples, the math won't work out right.

Potential overflow

When I used 10000 points and a grid size of 1000, the calculation overflowed because inbound exceeded the maximum integer value.  You might want to use a 64-bit value instead of an int.

pow(x,2)

I have a pet peeve against writing pow(x,2).  I think that x*x is both clearer and faster.  Check my profile and you'll see.

Rewrite

Here is a rewrite of your program with the following additional changes:


Sped up the inner loop by multiplying against a precomputed value in order to avoid costly divisions.
Cosmetic changes to spacing, braces, indentation, and location of variable declarations.
Added temp variables to simplify the lengthy equations.




#include &lt;stdio.h&gt;
#include &lt;stdlib.h&gt;
#include &lt;time.h&gt;
#include &lt;math.h&gt;
#include &lt;stdint.h&gt;

const double pi_exactly=3.14159265358979323846;

int main(void)
{
    int      m       = 0;
    int      g       = 0;
    uint64_t inbound = 0;

    srand(time(NULL));

    printf("Quantity of random points per square: \n");
    scanf("%d", &amp;m);

    printf("Gridsize: \n");
    scanf("%d", &amp;g);

    double rand_multiplier = 1.0 / ((double) RAND_MAX * g);

    for (int i=0; i&lt;g; i++) {
        double x  = (double) i / g;
        double x1 = (double) (i+1) / g;
        for (int j=0; j&lt;g; ++j) {
            double y  = (double) j / g;
            double y1 = (double) (j+1) / g;

            if (x1*x1 + y1*y1 &lt;= 1) {
                inbound += m;
                continue;
            } else if (x*x + y*y &gt; 1) {
                break;
            }

            for (int k=0; k&lt;m; k++) {
                double xr = x + rand() * rand_multiplier;
                double yr = y + rand() * rand_multiplier;

                if (xr*xr + yr*yr &lt;= 1)
                    inbound++;
            }
        }
    }


    double pi = (double) inbound / ((double) m * g*g)*4;
    double uncertainty = pi - pi_exactly;

    printf("The approximation of pi is: %lf\n", pi);
    printf("Uncertainty: %e\n", uncertainty);
    return 0;
}