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

typedef struct
{
  int   size;
  int   index;
  float sum;
  float *x;
} Cic;

typedef struct
{
  float recip;
  Cic *x;
  Cic *xx;
} Run_Stats;

typedef struct
{
  float a;
  float b;
  float y;
} Exp_Avg;

typedef struct
{
  float recip;
  Exp_Avg *x;
  Exp_Avg *xx;
} Exp_Stats;

Cic* cic_open(int size)
{
  Cic *s;

  if ((s = (Cic*)malloc(sizeof(Cic))) != NULL)
  {
    if ((s->x = (float*)calloc(size, sizeof(float))) == NULL)
    {
      free(s);
      s = 0;
    }
    else
    {
      s->size = size;
      s->index = 0;
      s->sum = 0.0;
    }
  }
  return s;
}

void cic_close(Cic *s)
{
  free(s->x);
  free(s);
}

float cic_eval(Cic *s, float x)
{
  int index;

  index = s->index;
  // subtract the oldest sample from the sum
  s->sum -= s->x[index];
  // overwrite the oldest sample with the newest one
  s->x[index++] = x;
  // add the newest sample to the sum
  s->sum += x;
  // update the index
  s->index = (index >= s->size) ? 0 : index;

  return s->sum;
}

Exp_Avg* expavg_open(int size)
{
  Exp_Avg *s;

  if ((s=malloc(sizeof(Exp_Avg))))
  {
    s->y = 0.0;
    s->a = size/(size+1.0);
    s->b = 1.0 - s->a;
  }
  return s;
}

void expavg_close(Exp_Avg *s)
{
  free(s);
}

float expavg_eval(Exp_Avg *s, float x)
{
  float y;

  y = s->a * s->y + s->b * x;
  s->y = y;
  return y;
}

Exp_Stats* expstat_open(int size)
{
  Exp_Stats *s = 0;

  if ((size > 3) && (s = (Exp_Stats*)malloc(sizeof(Exp_Stats))))
  {
    if(!(s->x = expavg_open(size)))
    {
      free(s);
      s = 0;
    }
    else if(!(s->xx = expavg_open(size)))
    {
      expavg_close(s->x);
      free(s);
      s = 0;
    }
    else
    {
      s->recip = 1.0/size;
    }
  }
  return s;
}


void expstat_close(Exp_Stats *s)
{
  expavg_close(s->x);
  expavg_close(s->xx);
  free(s);
}

void expstat(Exp_Stats *s, float x)
{
  expavg_eval(s->x, x);
  expavg_eval(s->xx, x*x);
}


float exp_mean(Exp_Stats *s)
{
  return s->x->y;
}

float exp_var(Exp_Stats *s)
{
  return s->xx->y - (s->x->y * s->x->y);
}

float exp_std(Exp_Stats *s)
{
  return sqrt(exp_var(s));
}

Run_Stats *runstat_open(int size)
{
  Run_Stats *s = 0;

  if ((size > 3) && (s = (Run_Stats*)malloc(sizeof(Run_Stats))))
  {
    if(!(s->x = cic_open(size)))
    {
      free(s);
      s = 0;
    }
    else if(!(s->xx = cic_open(size)))
    {
      cic_close(s->x);
      free(s);
      s = 0;
    }
    else
    {
      s->recip = 1.0/size;
    }
  }
  return s;
}

void runstat_close(Run_Stats *s)
{
  cic_close(s->x);
  cic_close(s->xx);
  free(s);
}

void runstat(Run_Stats *s, float x)
{
  cic_eval(s->x, x);
  cic_eval(s->xx, x*x);
}

float run_mean(Run_Stats *s)
{
  return s->x->sum * s->recip;
}

float run_var(Run_Stats *s)
{
  float mean;
  mean = run_mean(s);
  return (s->recip * s->xx->sum) - (mean*mean);
}

float run_std(Run_Stats *s)
{
  return sqrt(run_var(s));
}

void usage(char *prog)
{
  printf("Usage: %s <size> [<filename>]\n", prog);
  exit(-1);
}

int main(int argc, char **argv)
{
  Run_Stats *r;
  Exp_Stats *e;
  int size;
  FILE *infile = stdin;
  float x;
  int i;

  if ((argc < 2) || (argc > 3))
  {
    usage(argv[0]);
  }
  size = atoi(argv[1]);
  if (size < 3)
  {
    printf("The size argument must be greater than 2\n");
    return -1;
  }
  if ((argc == 3) && ((infile = fopen(argv[2], "r")) == NULL))
  {
    printf("Could not open %s for reading\n", argv[2]);
    return -1;
  }

  r = runstat_open(size);
  e = expstat_open(size);
  i = 0;
  while (fscanf(infile, "%f", &x) != EOF)
  {
    runstat(r,x);
    expstat(e,x);
  }
  printf("run avg=%f, run var=%f run std=%f\n", 
         run_mean(r), run_var(r), run_std(r));
  printf("exp avg=%f, exp var=%f exp std=%f\n", 
         exp_mean(e), exp_var(e), exp_std(e));
  runstat_close(r);
  expstat_close(e);
  return 0;
}