implicit none

integer, parameter :: num=1000
real :: x(num), y(num)
integer :: i, j, n
real :: xm, ym, sigx, sigy, cor
real :: a, b

  i=1
  do while( 1==1 )
    read(*,*, END=900) x(i), y(i)
    i=i+1
    if(i>num) then
      write(0,*) 'too many data. Change num'
      stop
    end if
  end do
900 n=i-1

  do i=1, n
   write(0,*) i, x(i), y(i)
  end do

 call correlation(x, y, n, xm, ym, sigx, sigy, cor)

  write(*,*) 'Mean x = ', xm, 'Sig x =', sigx
  write(*,*) 'Mean y = ', ym, 'Sig y =', sigy
  write(*,*) 'Correlation = ', cor

 call calc_ab( x, y, n, a, b)
  write(*,*) 'y=ax+b', 'a=', a, 'b=', b

end


subroutine correlation(x, y, n, xm, ym, sigx, sigy, cor)

 implicit none
 integer :: n, i
 real    :: x(n), y(n)
 real :: xm, ym, sigx, sigy, cor

 real(kind=8) :: sx, sy, sxy, sx2, sy2
  
 sx=0; sy=0; sxy=0; sx2=0; sy2=0


if( n > 2 ) then
 do i=1, n
  sx = sx +x(i)
  sy = sy +y(i)
 end do
  xm = sx/n; ym = sy/n
 do i=1, n
  sx2 = sx2 + (x(i)-xm)*(x(i)-xm)
  sy2 = sy2 + (y(i)-ym)*(y(i)-ym)
  sxy = sxy + (x(i)-xm)*(y(i)-ym)
 end do

  sigx = sqrt(sx2/n); sigy = sqrt(sy2/n)
  if( sigx > 0 .AND. sigy > 0 ) then
  cor = (sxy/n)/(sigx*sigy)
  end if
end if
end


subroutine calc_ab( x, y, n, a, b)
!   y = ax+b 
 implicit none
 integer :: n, i
 real :: x(n), y(n), a, b
 double precision :: sx2, sx, sxy, sy, d

 sx2 = 0 ;  sx = 0 ;  sxy = 0;  sy = 0

 do i=1,n
   sx2 = sx2 + x(i)*x(i) ! sx2 = sx2 + i*i
   sx = sx + x(i)        ! sx  = sx + i
   sxy = sxy + x(i)*y(i) !sxy = sxy + i*x(i)
   sy = sy + y(i)        !sy  = sy  + x(i)
 end do

 d = sx2 * n - sx * sx
 a = (n*sxy-sx*sy)/d
 b = (-sx*sxy+sx2*sy)/d
end







