! sigma_error.f90  show difference in errors

function stddeva(n, x) result(sdev)
  implicit none
  integer :: n
  real, dimension(*) :: x
  real :: sdev
  real :: sumx2, sumx, rn
  integer i

  rn = float(n)
  sumx2 = 0.0
  sumx = 0.0
  do i=1,n
    sumx = sumx + x(i)
    sumx2 = sumx2 + x(i)*x(i)
  end do
  sdev = sqrt((sumx2 - sumx*sumx/rn)/rn)
end function stddeva

function stddevb(n, x) result(sdev)
  implicit none
  integer :: n
  real, dimension(*) :: x
  real :: sdev
  real :: mean, sum2, rn
  integer i

  rn = float(n)
  mean = 0.0
  sum2 = 0.0
  do i=1,n
    mean = mean + x(i)
  end do
  mean = mean/rn
  do i=1,n
    sum2 = sum2 +(x(i)-mean)**2
  end do
  sdev = sqrt(sum2/rn)
end function stddevb


program sigma_error
  real :: sigma1, sigma2
  double precision :: sigma
  real, dimension(100) :: x
  integer :: n, i
  n = 100
  print *, 'sigma_error.f90 running, n=', n

  do i=1,n
    x(i) = float(i)
  end do
  print *, 'x 1 .. 100'

  sigma1 = stddeva(n, x)
  print *, 'sum x^2 - (sum x)^2/n  sigma1 = ', sigma1

  sigma2 = stddevb(n, x)
  print *, 'sum(x-mean)2/n  sigma2 = ', sigma2

  do i=1,n
    x(i) = 10.0*float(i)
  end do
  print *, 'x 10 .. 1000'

  sigma1 = stddeva(n, x)
  print *, 'sum x^2 - (sum x)^2/n  sigma1 = ', sigma1

  sigma2 = stddevb(n, x)
  print *, 'sum(x-mean)2/n  sigma2 = ', sigma2

  do i=1,n
    x(i) = 100.0*float(i)
  end do
  print *, 'x 100 .. 10000'

  sigma1 = stddeva(n, x)
  print *, 'sum x^2 - (sum x)^2/n  sigma1 = ', sigma1

  sigma2 = stddevb(n, x)
  print *, 'sum(x-mean)2/n  sigma2 = ', sigma2

  do i=1,n
    x(i) = 1000.0*float(i)
  end do
  print *, 'x 1000 .. 100000'

  sigma1 = stddeva(n, x)
  print *, 'sum x^2 - (sum x)^2/n  sigma1 = ', sigma1

  sigma2 = stddevb(n, x)
  print *, 'sum(x-mean)2/n  sigma2 = ', sigma2

  do i=1,n
    x(i) = 10000.0*float(i)
  end do
  print *, 'x 10000 .. 1000000'

  sigma1 = stddeva(n, x)
  print *, 'sum x^2 - (sum x)^2/n  sigma1 = ', sigma1

  sigma2 = stddevb(n, x)
  print *, 'sum(x-mean)2/n  sigma2 = ', sigma2

  do i=1,n
    x(i) = 10000.0+float(i)
  end do
  print *, 'x 10001 .. 10100'

  sigma1 = stddeva(n, x)
  print *, 'sum x^2 - (sum x)^2/n  sigma1 = ', sigma1

  sigma2 = stddevb(n, x)
  print *, 'sum(x-mean)2/n  sigma2 = ', sigma2

end program sigma_error