#include <stdio.h>
#include <math.h>

typedef unsigned long h32; // 32 bit hash
typedef unsigned long u32; // 32 bit unsigned int

#define yk_kMersenne   ((long) 0x7FFFFFFF) /* mersenne prime */
#define yk_kMultiplier ((long) 48271)      /* prime multiplier */
#define yk_kQuotient   ((long) 44488)      /* mersenne / multiplier */
#define yk_kRemainder  ((long) 3399)       /* mersenne % multiplier */

h32 rand_j(h32 inSelf) // standard Park and Miller pseudo rand number
{
  if ( inSelf ) { // is self non-zero (we will not return zero)?
    long high = (long) inSelf / yk_kQuotient;
    long low = (long) inSelf % yk_kQuotient;
    long test = (yk_kMultiplier * low) - (yk_kRemainder * high);
    inSelf = (h32) (( test > 0 )? test : test + yk_kMersenne);
  }
  else
    inSelf = 1; // never return zero
    
  return inSelf;
}

u32 permute(void** v, int length, u32 seed)
{
  for (int i = length; i > 0; i--) {
    seed = rand_j(seed); // next random number
    int pick = (int) (seed % (i+1)); // random index in range 0..i
    if (pick != i) { // actually need to swap
        void* tmp = v[i]; v[i] = v[pick]; v[pick] = tmp;
    }
  }
  return seed;
}

double // random float in exclusive interval (0,1)
double_j(h32 inSelf, h32* outSeed)
{
  if ( inSelf ) { /* is self non-zero (we will not return zero)? */
    long high = (long) inSelf / yk_kQuotient;
    long low = (long) inSelf % yk_kQuotient;
    long test = (yk_kMultiplier * low) - (yk_kRemainder * high);
    inSelf = (h32) (( test > 0 )? test : test + yk_kMersenne);
  }
  else
    inSelf = 1; /* never return zero */
    
  if ( outSeed )
    *outSeed = inSelf;
  return ((double) inSelf) / ((double) yk_kMersenne);
}

class Statistic { // variance for single statistic
public:
    long long  m_len; // number of points observed
    double     m_sum; // sum of all points
    double     m_squares; // sum of squared points
public:
    Statistic() : m_len(0), m_sum(0.0), m_squares(0.0) { }
public:
    void add(u32 x) {
        ++m_len;
        double real = x;
        m_sum += real;
        m_squares += real * real;
    }
    void add(double x) {
        ++m_len;
        m_sum += x;
        m_squares += x * x;
    }
    double mean() const { return (m_len>1)? m_sum/m_len : m_sum; }
    double variance() const { // square of standard deviation
        if (m_len > 1) {
            double mu = this->mean();
            double n_mu_sq = m_len * (mu * mu);
            return (m_squares - n_mu_sq)/(m_len - 1);
        }
        return 0.0;
    }
};

class Normal { // normal distribution
protected:
    h32     m_seed; // random number generator seed
    double  m_mean; // average value generated
    double  m_sdev; // standard deviation
    double  m_now;  // last value returned by operator*()
    
    double  m_y2;   // second value generated by next()
    int     m_use_last; // whether to use m_y2
public:
    Normal(u32 seed, double mean, double sdev);
    h32 seed() const { return m_seed; }
    
    double next();
    double operator++() { return (m_now = next()); }
    double operator*() { return m_now; }
};

Normal::Normal(u32 seed, double mean, double sdev)
    : m_seed(seed)
      , m_mean(mean)
      , m_sdev(sdev)
      , m_now(mean)
      , m_y2(0.0)
      , m_use_last(0)
{
}

double Normal::next()
{
  double x1, x2, w, y1;
  //static float y2;
  //static int use_last = 0;

  if (m_use_last)
  {
    y1 = m_y2;
    m_use_last = 0;
  }
  else
  {
    do {
      x1 = 2.0 * double_j(m_seed, &m_seed) - 1.0;
      x2 = 2.0 * double_j(m_seed, &m_seed) - 1.0;
      w = x1 * x1 + x2 * x2;
    } while ( w >= 1.0 );

    w = ::sqrt( (-2.0 * ::log( w ) ) / w );
    y1 = x1 * w;
    m_y2 = x2 * w;
    m_use_last = 1;
  }

  return m_mean + (y1 * m_sdev);
}

static void stars(int n)
{
    int top = n / 4;
    for (int i = 0; i < top; i++)
        fputc('#', stdout);

    int last = n % 4;
    if (last) {
        static const char* odd = " -+*"; // 0 1 2 3
        int c = odd[last];
        fputc(c, stdout);
    }
}

#define TMP_ARRY_LEN 2048
#define TMP_HILL_LEN 30

int variance_examples()
{
    int hill[ TMP_HILL_LEN + 1 ];
    u32 k;
    for (k = 0; k < TMP_HILL_LEN; k++)
        hill[k] = 0;

    double ary[ TMP_ARRY_LEN + 1 ];

    Statistic st;
    Normal nm(31555, 1500.0, 450.0);
   
    for (int i = 0; i < TMP_ARRY_LEN; i++) {
        double point = ++nm;
        ary[ i ] = point;
        st.add(point);
        u32 j = (u32) (point / 100.0);
        if (j < TMP_HILL_LEN)
            ++hill[j];
    }
    double avg = st.mean();
    double var = st.variance();
    double sd = (var<0)? ::sqrt(-var): ::sqrt(var);

    fprintf(stdout, "avg='%lf' var='%lf' sd='%lf' cv='%lf'\n", avg, var, sd, sd/avg);

    for (k = 0; k < TMP_HILL_LEN; k++) {
        fprintf(stdout, "\n%04d-%04d: ", k*100, (k+1)*100);
        stars(hill[k]);
    }

    fputc('\n', stdout);
    fflush(stdout);

    return 0;
}

int main (int argc, char * const argv[]) {
    variance_examples();
    return 0;
}

/*
[Session started at 2005-11-22 07:34:07 -0800.]
avg='1502.054369' var='200076.060876' sd='447.298626' cv='0.297791'

0000-0100: +
0100-0200: *
0200-0300: *
0300-0400: #*
0400-0500: ###-
0500-0600: ###+
0600-0700: #######-
0700-0800: ###########-
0800-0900: #############
0900-1000: ##########################*
1000-1100: ##########################
1100-1200: #####################################*
1200-1300: ###################################*
1300-1400: ##########################################
1400-1500: ############################################-
1500-1600: ##########################################-
1600-1700: ##############################################*
1700-1800: #########################################*
1800-1900: ###############################*
1900-2000: ########################
2000-2100: ###########################
2100-2200: #################
2200-2300: #######*
2300-2400: #######*
2400-2500: ####
2500-2600: ##+
2600-2700: ##+
2700-2800: #-
2800-2900: +
2900-3000: +

boxmuller has exited with status 0.
*/