c*******************************************
      program least_squares
c*******************************************
      implicit double precision(a-h,o-z)
      parameter (ndata=100,npol=10)
      dimension x(ndata),y(ndata)
      dimension a(npol,npol),b(npol),coef(npol)
c
      open(1,file='data')
      do i=1,ndata
        read(1,*,end=10)x(i),y(i)
      enddo
      stop     'Increase NDATA'
  10  npoints=i-1
      close(1)
c
      write(*,*)' Enter polynomial degree '
      read(*,*)mpol
      if(mpol.ge.npol) stop 'Increase NPOL'
c
      do i=1,mpol+1
      do j=i,mpol+1
        call fill_lhs(ndata,npoints,x,i,j,res)
        a(i,j)=res
        a(j,i)=res
      enddo
        call fill_rhs(ndata,npoints,x,y,i,b(i))
      enddo
c
      call gauss(npol,mpol+1,a,b,coef)
c
      open(1,file='approx')
      do i=1,npoints
        app=0.d0
      do j=1,mpol+1
        app=app+coef(j)*x(i)**(j-1)
      enddo
      error=100.*(app-y(i))/y(i)
      write(1,'(2x,f7.0,2x,f12.5,2x,f10.5)')x(i),app,error
      enddo
      close(1)
c
      end
c
c
c***********************************************************************
      subroutine fill_lhs(ndata,npoints,x,i,j,res)
c***********************************************************************
      implicit double precision(a-h,o-z)
      dimension x(ndata)
      res=0.d0
      expo=dfloat(i+j-2)
      do m=1,npoints
        res = res + x(m)**expo
      enddo
      return
      end
c
c***********************************************************************
      subroutine fill_rhs(ndata,npoints,x,y,i,res)
c***********************************************************************
      implicit double precision(a-h,o-z)
      dimension x(ndata),y(ndata)
      res=0.d0
      do m=1,npoints
        res = res + y(m)*x(m)**dfloat(i-1)
      enddo
      return
      end
c
c
c***********************************************************************
      subroutine gauss(kdim,n,a,b,x)
c***********************************************************************
      implicit double precision (a-h,o-z)
      dimension a(kdim,kdim),b(kdim),x(kdim)
      eps=1.d-8
      do k=1,n-1
      if (dabs(a(k,k)).lt.eps)   stop   'Divide by Zero !!!!'
      do i=k+1,n 
      factor=a(i,k)/a(k,k)
      do j=k+1,n
      a(i,j)=a(i,j)-factor*a(k,j)
      enddo
      b(i)  =b(i)  -factor*b(k)
      enddo
      enddo
c
      x(n)=b(n)/a(n,n)
      do i=n-1,1,-1
      sum=0.d0
      do j=i+1,n
      sum=sum+a(i,j)*x(j)
      enddo
      x(i)=(b(i)-sum)/a(i,i)
      enddo
c
      return
      end