-- sigma_error.adb  show error for different methods

with Ada.Text_IO; use Ada.Text_IO;
with Ada.Numerics.Elementary_Functions;
use  Ada.Numerics.Elementary_Functions;

procedure Sigma_Error is
  type Float_Arr is array(integer range <>) of float;
  sigma1, sigma2 : float;
  n : integer := 100;
  x : Float_Arr(1..n);

  function stddeva(n: integer; x : Float_Arr) return Float is
    sdev, sumx2, sumx, rn : float;
  begin
    rn := float(N);
    sumx2 := 0.0;
    sumx := 0.0;
    for i in 1..n loop
      sumx := sumx + x(i);
      sumx2 := sumx2 + x(i)*x(i);
    end loop;
    sdev := sqrt((sumx2 - sumx*sumx/rn)/rn);
    return sdev;
  end stddeva;

  function stddevb(n: integer; x :Float_Arr) return Float is
    sdev, mean, sum2, rn : float;
  begin
    rn := float(N);
    mean := 0.0;
    sum2 := 0.0;
    for i in 1..n loop
      mean := mean + x(i);
    end loop;
    mean := mean/rn;
    for i in 1..n loop
     sum2 := sum2 +(x(i)-mean)*(x(i)-mean);
    end loop;
    sdev := sqrt(sum2/rn);
    return sdev;
  end Stddevb;

begin

  Put_Line("sigma_error.adb running, n="&integer'image(N));

  for i in 1..n loop
    x(i) := 1.0+float(I);
  end loop;
  Put_Line("x 1 .. 100 ");

  sigma1 := stddeva(n, x);
  Put_Line("sum x^2 - (sum x)^2/n  sigma1 := "&float'image(sigma1));

  sigma2 := stddevb(n, x);
  Put_Line("sum(x-mean)2/n  sigma2 := "&float'image(sigma2));

  for i in 1..n loop
    x(i) := 1.0+10.0*float(i);
  end loop;
  Put_line("x 10 .. 1000 ");

  sigma1 := stddeva(n, x);
  Put_Line("sum x^2 - (sum x)^2/n  sigma1 := "&float'image(sigma1));

  sigma2 := stddevb(n, x);
  Put_Line("sum(x-mean)2/n  sigma2 := "&float'image(sigma2));

  for i in 1..n loop
    x(i) := 1.0+100.0*float(i);
  end loop;
  Put_Line("x 100 .. 10000 ");

  sigma1 := stddeva(n, x);
  Put_Line("sum x^2 - (sum x)^2/n  sigma1 := "&float'image(sigma1));

  sigma2 := stddevb(n, x);
  Put_Line("sum(x-mean)2/n  sigma2 := "&float'image(sigma2));

  for i in 1..n loop
    x(i) := 1.0+1000.0*float(i);
  end loop;
  Put_Line("x 1000 .. 100000 ");

  sigma1 := stddeva(n, x);
  Put_Line("sum x^2 - (sum x)^2/n  sigma1 := "&float'image(sigma1));

  sigma2 := stddevb(n, x);
  Put_Line("sum(x-mean)2/n  sigma2 := "&float'image(sigma2));

  for i in 1..n loop
    x(i) := 1.0+10000.0*float(i);
  end loop;
  Put_Line("x 10000 .. 1000000 ");

  sigma1 := stddeva(n, x);
  Put_Line("sum x^2 - (sum x)^2/n  sigma1 := "&float'image(sigma1));

  sigma2 := stddevb(n, x);
  Put_Line("sum(x-mean)2/n  sigma2 := "&float'image(sigma2));

  for i in 1..n loop
    x(i) := 1.0+10000.0+float(i);
  end loop;
  Put_Line("x 10001 .. 10100 ");

  sigma1 := stddeva(n, x);
  Put_Line("sum x^2 - (sum x)^2/n  sigma1 := "&float'image(sigma1));

  sigma2 := stddevb(n, x);
  Put_Line("sum(x-mean)2/n  sigma2 := "&float'image(sigma2));

end Sigma_Error;