/* sigma_error.c show error for different methods */

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

float stddeva(int n, float x[])
{
  float sdev, sumx2, sumx, rn;
  int i;

  rn = (float)n;
  sumx2 = 0.0;
  sumx = 0.0;
  for(i=0; i<n; i++)
  {
    sumx = sumx + x[i];
    sumx2 = sumx2 + x[i]*x[i];
  }
  sdev = sqrt((sumx2 - sumx*sumx/rn)/rn);
  return sdev;
}

float stddevb(int n, float x[])
{
  float sdev, mean, sum2, rn;
  int i;

  rn = (float)n;
  mean = 0.0;
  sum2 = 0.0;
  for(i=0; i<n; i++)
    mean = mean + x[i];
  mean = mean/rn;
  for(i=0; i<n; i++)
    sum2 = sum2 +(x[i]-mean)*(x[i]-mean);
  sdev = sqrt(sum2/rn);
  return sdev;
}


int main(int argc, char *argv[])
{
  float sigma1, sigma2;
  float x[100];
  int n, i;
  n = 100;
  printf("sigma_error.c running, n=%d \n", n);

  for(i=0; i<n; i++)
    x[i] = (float)(i+1);
  printf("x 1 .. 100 \n");

  sigma1 = stddeva(n, x);
  printf("sum x^2 - (sum x)^2/n  sigma1 = %g \n", sigma1);

  sigma2 = stddevb(n, x);
  printf("sum(x-mean)^2/n  sigma2 = %g \n", sigma2);

  for(i=0; i<n; i++)
    x[i] = 10.0*(float)(i+1);
  printf("x 10 .. 1000 \n");

  sigma1 = stddeva(n, x);
  printf("sum x^2 - (sum x)^2/n  sigma1 = %g \n", sigma1);

  sigma2 = stddevb(n, x);
  printf("sum(x-mean)2/n  sigma2 = %g \n", sigma2);

  for(i=0; i<n; i++)
    x[i] = 100.0*(float)(i+1);
  printf("x 100 .. 10000 \n");

  sigma1 = stddeva(n, x);
  printf("sum x^2 - (sum x)^2/n  sigma1 = %g \n", sigma1);

  sigma2 = stddevb(n, x);
  printf("sum(x-mean)^2/n  sigma2 = %g \n", sigma2);

  for(i=0; i<n; i++)
    x[i] = 1000.0*(float)(i+1);
  printf("x 1000 .. 100000 \n");

  sigma1 = stddeva(n, x);
  printf("sum x^2 - (sum x)^2/n  sigma1 = %g \n", sigma1);

  sigma2 = stddevb(n, x);
  printf("sum(x-mean)^2/n  sigma2 = %g \n", sigma2);

  for(i=0; i<n; i++)
    x[i] = 10000.0*(float)(i+1);
  printf("x 10000 .. 1000000 \n");

  sigma1 = stddeva(n, x);
  printf("sum x^2 - (sum x)^2/n  sigma1 = %g \n", sigma1);

  sigma2 = stddevb(n, x);
  printf("sum(x-mean)^2/n  sigma2 = %g \n", sigma2);

  for(i=0; i<n; i++)
    x[i] = 10000.0+(float)(i+1);
  printf("x 10001 .. 10100 \n");

  sigma1 = stddeva(n, x);
  printf("sum x^2 - (sum x)^2/n  sigma1 = %g \n", sigma1);

  sigma2 = stddevb(n, x);
  printf("sum(x-mean)^2/n  sigma2 = %g \n", sigma2);

} /* end sigma_error.c */