! peval.f90 Horners method for evaluating a polynomial 

function peval(n, x, c) result (y)
  ! Horners method for evaluating a polynomial   
  ! an nth order polynomial has n+1 coefficients 
  ! stored in c(0) through c(n)                  
  ! y = c(0) + c(1)*x + c(2)*x^2 +...+ c(n)*x^n  
  implicit none
  integer, intent(in) :: n
  double precision, intent(in) :: x
  double precision, dimension(0:*), intent(in) :: c
  double precision :: y
  integer i

  if(n<=0) y = c(0); return
  y = c(n)*x
  do i=n-1, 1, -1
    y = (c(i)+y)*x
  end do
  y = y+c(0)
end function peval

program test_peval
  double precision :: y
  double precision, dimension(0:11) :: c, r
  integer :: i, j, n

  interface
    function peval(n, x, c) result (y)
      integer, intent(in) :: n
      double precision, intent(in) :: x
      double precision, dimension(0:*), intent(in) :: c
      double precision :: y
    end function peval
  end interface

  print *, 'peval.f90 running'
  ! test peval by generation roots of a polynomial, 
  ! then generating the polynomial, then evaluating 
  do n=3,10  ! n is degree of polynomial 
    print *, ' '
    print *, 'y = c(0) + c(1)*x + c(2)*x^2 +...+ c(n)*x^n, for n', n
    print *, 'roots are:'
    r(0) = 0.1D0 
    print *, 'r(0)=', r(0) 
    r(1) = 2.5D0
    print *, 'r(1)=', r(1) 
    do i=2,n-1
      r(i) = r(i-1)+r(i-2) 
      print *, 'r(',i,')=', r(i) 
    end do
    ! generate polynomial coefficients 
    c(0) = r(0)*r(1) 
    c(1) = -(r(0)+r(1)) 
    c(2) = 1.0 
    do j=2,n-1
      do i=j+1,1,-1
        c(i) = c(i-1)
      end do 
      c(0) = 0.0 
      do i=0,j 
        c(i) = c(i) - r(j)*c(i+1)
      end do 
    end do
    print *, 'polynomial coefficients are:' 
    do i=0,n
      print *, 'c(',i,')=', c(i)
    end do 
    print *, 'polynomial evaluated at each root:'
    do i=0, n-1
      y = peval(n, r(i), c) 
      print *, 'r(',i,')=',r(i),', peval(',n,',r(',i,'),c)=',y
    end do
  end do
end program test_peval