#include <stdio.h>
#include <stdlib.h>
#include <sys/time.h>
#include <math.h>
#define PTS 10000000
#define MXSTRIDE 20
void main() {
  unsigned long i, k, istride, msec;
  double *data, cnorm, mean, stddev, u1, u2, xk, deltmean, recipk, rate;
  double two_pi = 2.0 * M_PI;
  struct timeval tm1, tm2;
  data = malloc(PTS*MXSTRIDE*sizeof(double));
  cnorm = 1.0 / (RAND_MAX + 1.0);
  for (i = 0; i <= PTS*MXSTRIDE; i++) {
    data[i] = (rand() + 1) * cnorm;
  }
  for (istride = 1; istride <= MXSTRIDE; istride++) {
    k = 1;
    stddev = 0;
    mean = 0.;
    u1 = (rand() + 1) * cnorm;  // need one extra data point
    gettimeofday(&tm1, NULL);
    for (i = 0; i < PTS*istride; i += istride) {
      // convert data to normal distribution using Box-Muller
      u2 = u1;
      u1 = data[i];
      xk = sqrt(-2.0 * log(u1)) * cos(two_pi * u2);
      // compute the cumulative mean and standard deviation*
      deltmean = xk - mean;
      recipk = 1.0/k;
      mean = mean + deltmean*recipk;
      stddev = stddev + (k - 1)*(deltmean*deltmean)*recipk;
      k += 1;
    }
    gettimeofday(&tm2, NULL);
    stddev = stddev / k;
    msec = 1000 * (tm2.tv_sec - tm1.tv_sec) + \
        (tm2.tv_usec - tm1.tv_usec) / 1000;
    rate = sizeof(double)*PTS*(1000.0/msec) / (1024*1024);
    printf("stride %2lu, mean %9.6f, stddev %8.6f ", istride, mean, stddev);
    printf("time %3lu msec, rate %4.0f MB/s\n", msec, rate);
  }
// *https://citeseerx.ist.psu.edu/document?repid=rep1&type=pdf&doi=d28cb3d22420bc144898b71373843f9114139206
// Formulas originally found at http://www.cs.berkeley.edu/~mhoemmen/cs194/Tutorials/variance.pdf
}