#include <iostream>
#include <cmath>
using namespace std;

int median( int, int, int );

/*
 *  This program demonstrates that chosing the median-of-three
 *  divides the interval into sub-intervals of approximately
 *  width 5/16 and 11/16.
 *
 *             1/2
 *            /
 *           |                       5
 *        2  |     6 x (1 - x) dx = -- = 0.3125
 *           |                      16
 *          /
 *            0
 *
 *  This program estimates the width of the smaller 
 *  interval by randomly chosing one billion evenly
 *  distributed integers on [0, 2^20 - 1] and finding
 *  the width of the smaller interval.
 *
 *  After 394 s, this program prints the value
 *      Mean:                           0.31249974682817804883
 *      Standard Deviation:             0.1218350096311656694
 *      Standard Deviation of the mean: 1.2183500963116567517e-06
 *
 *  The good news is that the average is over 2.5 standard
 *  deviations from pathalogical case:  splitting the
 *  array into subintervals of size 1 and N - 2.
 */

int const N = 1048575;         // 2^20 - 1 
long const M = 10000000000;    // 10 billion

int main() {
	srand48( 2009 );
	cout.precision( 20 );

	double ratio = 0.0;
	double variance = 0.0;

	for ( long i = 0; i < M; ++i ) {
		int med = median( lrand48() & N, lrand48() & N, lrand48() & N );
		double dmed = static_cast<double>( min( med, N - med ) );
		ratio += dmed;
		variance += dmed*dmed;
	}

	double mean = ratio/N/M;

	// Print out the average
	cout << "Mean:                           "
	     << mean << endl;
	cout << "Standard Deviation:             "
	     << sqrt( variance/N/N/M - mean*mean ) << endl;
	cout << "Standard Deviation of the mean: "
	     << sqrt( variance/N/N/M/M - mean*mean/M ) << endl;

	return 0;
}

/*
 * median( a, b, c )
 *
 *  Return the median of the three arguments.
 *    min( a, b, c ) <= median( a, b, c ) <= max( a, b, c )
 */

int median( int a, int b, int c ) {
	return b < c ?
		( a < b ? b : ( a < c ? a : c ) ) :
		( a < c ? c : ( a < b ? a : b ) );
}