//  /*--------------------------------------------------------*\
//  |   Main Program: Monte Carlo estimate of pi
//  |*--------------------------------------------------------*|
//  |   Date:         14 August 1996 (revised 24 November 2006)
//  |   Author:       Bruce W. Weide
//  |   
//  |   Brief User's Manual:
//  |   Asks for a positive integer, which is the number of
//  |   points to use for a Monte Carlo estimate of pi, then
//  |   reports the resulting estimate of pi, and then quits.
//  |   
//  \*--------------------------------------------------------*/

///-------------------------------------------------------------
/// Global Context ---------------------------------------------
///-------------------------------------------------------------

#include "RESOLVE_Foundation.h"

#include "CI/Random/Uniform_Generator_1.h"

///-------------------------------------------------------------
/// Interface --------------------------------------------------
///-------------------------------------------------------------

program_body main ()
{
    object Character_IStream input;
    object Character_OStream output;
    object Random_Uniform_Generator_1 rn;
    object Integer n, pts_in_U, pts_in_Q;
    object Real pi_estimate;

    // Open input and output streams
    
    input.Open_External ("");
    output.Open_External ("");

    // Ask user for number of points
    
    output << "Number of points: ";
    input >> n;

    // Create counts for estimate

    while (pts_in_U < n)
    /*!
	preserves n
	alters rn, pts_in_U, pts_in_Q
	maintains
	    0 <= pts_in_Q <= pts_in_U <= n  and
	    pts_in_Q = [number of points that have been
	                in the quarter-circle]  and
	    pts_in_U = [number of points that have been
	                in the unit square]
	decreases
	    n - pts_in_U
    !*/
    {
        object Real x, y;

        // Generate a random point in the unit square U
    
        rn.Generate_Next ();
        x = rn.Uniform_Real (0.0, 1.0);
        rn.Generate_Next ();
        y = rn.Uniform_Real (0.0, 1.0);
        pts_in_U++;

        // Check if point is also in the quarter circle Q
    
        if (x * x + y * y < 1.0)
        {
            pts_in_Q++;
        }
    }

    // Estimate pi
    
    pi_estimate = 4.0 * To_Real (pts_in_Q) / To_Real (n);
    output << "pi is about equal to " << pi_estimate << "\n";

    // Close input and output streams
    
    input.Close_External ();
    output.Close_External ();
}