#include <iostream>
#include <utility>
#include <algorithm>
#include <ctime>
#include <cstdlib>
#include <vector>

typedef std::vector< double > DataVec;

double drnd( void )
{
    double aa = (double)::random();
    return aa / (double) 0x7FFFFFFFU;
}

double computeMean( DataVec& data ) {

    double mean = 0;

    for ( int i = 0; i < data.size(); i++ ) {
        mean += data[ i ];
    }

    return mean / data.size();
}

double classicalVariance( DataVec& data ) {

    double mean = computeMean( data );
    double variance = 0;

    for ( int i = 0; i < data.size(); i++ ) {
        variance += ( data[ i ] - mean ) * ( data[ i ] - mean );
    }

    return variance / data.size();

}

double provisionalMeans( DataVec& data ) {

    assert( data.size() );

    // compute first step, so we can use our first-order difference eqns
    double preMean = data[ 0 ];
    double preVar = 0.0;

    // init our running counts to pre values in case we have only one element
    double curMean = preMean;
    double curVar = preVar;

    // compute a running update of the mean and variance
    for ( int i = 1; i < data.size(); i++ ) {
        curMean = preMean + ( ( data[ i ] - preMean ) / ( i + 1 ) );
        curVar = preVar + ( ( data[ i ] - preMean ) * ( data[ i ] - curMean ) );
        preMean = curMean;
        preVar = curVar;
    }

    // kazaam!
    return curVar / data.size();

}

int main( int argc, char* argv[] ) {

    // our data vector
    DataVec data;

    // build "random" input data
    srand48( time( 0 ) );
    srand( time( 0 ) );
    static const int NUM_ELEMENTS = 800000;
    data.reserve( NUM_ELEMENTS );
    for ( int i = 0; i < NUM_ELEMENTS; i++ ) {
        data.push_back( drnd() );
    }

    std::cout << "Classical:   " << classicalVariance( data ) << std::endl;
    std::cout << "Provisional: " << provisionalMeans( data ) << std::endl;

}