Estimating Pi with the Monte Carlo method, loop appears to be terminating early

Question

I am attempting to write a program that estimates Pi based on the Monte Carlo method via a random number generator. I am attempting to estimate Pi within 1, 2, 3, 4, 5, and 6 digits of accuracy and have the program print to the screen how many points it took to get within .1 digit of Pi, then .01 digits of Pi and so forth all the way until .000001 digits of Pi. I am allowing the user to input an amount of trials they would like to run, so it will print "Trial 1, 2, 3, 4" etc. with all of the information I listed above. I am stuck on one last bit, and that is having it loop back through the calculations (it will not print more than just trial 1). Though I am not getting a message that the program has terminated, so I can not tell if it is my while loop failing or my nested for loops. Please help! :)

I have attempted switching around the for loops as well as trying different varying if statements. This is the closest I have gotten it to running the way I would like with exception of allowing the user to run multiple trials.

#include "pch.h"
#include <iostream> //need this by default for cin
#include <math.h> //includes math functions
#include <cmath> //includes basic math 
#include <cfloat> //includes floating point numbers
#include <iomanip> //includes setprecision for decimal places
#include <cstdlib> //needed for rand and srand functions
#include <ctime> //needed for time function used to seed generator
using namespace std;

int main()
{
cout << "The purpose of this program is to estimate pi using the monte 
carlo method and a random number generator" << endl << endl;

unsigned seed = time(0);
srand(seed);

float radius;
int trialcount = 0;
int trials;
float accuracy;
const float pi = 3.14159265;
float randpi = 0;
int squarecount = 0;
int circlecount = 0;
float x;
float y;
int n;


cout << "The value of PI can be found as the ratio of areas of a circle of radius r located within a square of side 2r" << endl;
cout << "This program runs a MonteCarlo Simulation that generates numbers located randomly within a square" << endl;
cout << "The count of values within the square and the count of numbers within the circle approximate their areas" << endl;
cout << "An input value of radius determines the size of the circle and square" << endl;
cout << "The user specifies how many trials or test runs are desired" << endl << endl;

cout << "The true value of PI to 8 decimal places is 3.14159265" << endl << endl;

cout << "Input a value for radius: "; 
cin >> radius;
cout << endl;
cout << "How many trials would you like? ";
cin >> trials;
cout << endl << endl;

cout << "Square count gives the Total number of random samples (they are within the square)" << endl;
cout << "Circle count gives the number of random samples that also fall within the circle" << endl << endl;


while (trialcount != trials)
{
    accuracy = .1;
    cout << "Trial " << trialcount + 1 << endl;
    cout << "Accuracy \t\t" << "Square Count \t\t" << "Circle Count \t\t" << "Pi" << endl << endl;

    for (int j = 0; randpi >= pi - accuracy || randpi <= pi + accuracy; j++)
    {
        cout << setprecision(6) << fixed << accuracy << " \t\t" << squarecount << " \t\t" << circlecount << " \t\t" << randpi << endl << endl;
        accuracy = accuracy / 10;

        for (int i = 0; randpi >= pi + accuracy || randpi <= pi - accuracy; i++)
        {
            x = (float)(rand());
            x = (x / 32767) * radius;
            y = (float)(rand());
            y = (y / 32767) * radius;

            squarecount++;

            if ((x * x) + (y * y) <= (radius * radius))
            {
                circlecount++;
            }

            randpi = float(4 * circlecount) / squarecount;

        }
    }

    trialcount++;

}


}

Answer 1

Problems I see:

Problem 1

The first for loop does not make any sense. If you want to make sure that you use accuracies of 0.1, 0.01, 0.001, etc. you just need a simple for loop. The following should do:

for ( int j = 0; j < 6; ++j )
{
    ...
}

Problem 2

x and y values are computed incorrectly. You want to make sure that their values are less than or equal to radius. However, when you use:

x = (float)(rand());
x = (x / 32767) * radius;
y = (float)(rand());
y = (y / 32767) * radius;

they are not guaranteed to be less than or equal to radius. They will be more than radius more often than they will not. You need to use

x = (float)(rand() % 32768);
x = (x / 32767) * radius;
y = (float)(rand() % 32768);
y = (y / 32767) * radius;

Problem 3

You need to reset the values of randpi, squarecount, and circlecount in every iteration of the inner for loop. Otherwise, your computations will be affected by computations from the previous iteration.

The outer for loop must start with:

for (int j = 0; j < 6; j++)
{
   accuracy /= 10;
   randpi = 0;
   squarecount = 0;
   circlecount = 0;

Problem 4

The inner for loop must be constrained to only run upto a certain number of times. If for some reason the accuracy is not achieved, you want to make sure you don't overflow i. For example:

int stopAt = (INT_MAX >> 8);
for (int i = 0; (randpi >= pi + accuracy || randpi <= pi - accuracy) && i < stopAt; i++)

For machines that use 32 bit ints, which is the most common in practice today, you won't run the loop any more than 0x7FFFFF (8388607 in decimal) times.

This is the core problem in your code. Your computations don't converge some times and you don't make sure you exit after a certain number of iterations of the loop.

Further improvement

You don't need radius as a variable in your program. You can compute x and y as:

x = (float)(rand() % 32768);
x = (x / 32767);
y = (float)(rand() % 32768);
y = (y / 32767);

and change the logic to check whether this is a point within the circle to

if ((x * x) + (y * y) <= 1.0 )

You should also attempt to define variables only in the scopes where you need them. This will make sure that you don't end up using stale values from a previous run of the iteration.

Revised program

The following revised program works for me.

#include <iostream> //need this by default for cin
#include <math.h> //includes math functions
#include <cmath> //includes basic math 
#include <cfloat> //includes floating point numbers
#include <iomanip> //includes setprecision for decimal places
#include <cstdlib> //needed for rand and srand functions
#include <ctime> //needed for time function used to seed generator
#include <climits> 

using namespace std;

int main()
{
   cout << "The purpose of this program is to estimate pi using the monte "
      "carlo method and a random number generator" << endl << endl;

   unsigned seed = time(0);
   srand(seed);

   int trialcount = 0;
   int trials;
   float accuracy;
   const float pi = 3.14159265;


   cout << "The value of PI can be found as the ratio of areas of a circle of radius r located within a square of side 2r" << endl;
   cout << "This program runs a MonteCarlo Simulation that generates numbers located randomly within a square" << endl;
   cout << "The count of values within the square and the count of numbers within the circle approximate their areas" << endl;
   cout << "An input value of radius determines the size of the circle and square" << endl;
   cout << "The user specifies how many trials or test runs are desired" << endl << endl;

   cout << "The true value of PI to 8 decimal places is 3.14159265" << endl << endl;

   cout << endl;
   cout << "How many trials would you like? ";
   cin >> trials;
   cout << endl << endl;

   cout << "Square count gives the Total number of random samples (they are within the square)" << endl;
   cout << "Circle count gives the number of random samples that also fall within the circle" << endl << endl;


   while (trialcount != trials)
   {
      accuracy = 0.1;
      cout << "Trial " << trialcount + 1 << endl;
      cout << "Accuracy \t\t" << "Square Count \t\t" << "Circle Count \t\t" << "Pi" << endl << endl;

      for (int j = 0; j < 6; j++)
      {
         accuracy /= 10;
         float randpi = 0;
         int squarecount = 0;
         int circlecount = 0;

         int stopAt = (INT_MAX >> 8);
         for (int i = 0; (randpi >= pi + accuracy || randpi <= pi - accuracy) && i < stopAt; i++)
         {
            float x = ((float)(rand() % 32768) / 32767);
            float y = ((float)(rand() % 32768) / 32767);

            squarecount++;

            if ((x * x) + (y * y) <= 1.0 )
            {
               circlecount++;
            }

            randpi = float(4 * circlecount) / squarecount;
         }

         cout << setprecision(8) << fixed << accuracy << " \t\t" << squarecount << " \t\t" << circlecount << " \t\t" << randpi << endl << endl;
      }

      trialcount++;
   }
}

See it working at https://ideone.com/laF27X.