program main
implicit none
real*8,allocatable, dimension(:) :: dati_x, dati_y
integer :: n_dati, i
real*8 :: pn, x

open(10,file='dati_lagrange.dat',status='OLD')
read(10,*) n_dati
n_dati = n_dati -1 
allocate(dati_x(0:n_dati),dati_y(0:n_dati))
do i=0,n_dati
  read(10,*) dati_x(i), dati_y(i)
enddo
close(10)

write(*,*)'Inserisci il nodo x dove calcolare Pn(x)'
write(*,'(1x,$,a)')'x='
read(*,*)x
write(*,*)'Pn(x)=',pn(x,dati_x,dati_y,n_dati)
deallocate(dati_x,dati_y)
end program main

function pn(x,dati_x,dati_y,n)
implicit none
real*8 :: pn, x, lagrange, sum
integer :: n, i
real*8, dimension(0:n) :: dati_x, dati_y
sum = 0.d0
do i=0,n
  sum = sum + dati_y(i)*lagrange(i,x,dati_x,n)
enddo
pn = sum
return
end function pn

function lagrange(i,x,dati_x,n)
implicit none
real*8 :: lagrange, x, prod
integer :: i, n, j
real*8, dimension(0:n) :: dati_x
prod = 1.d0
do j=0,n
  if (i /= j) then
    prod = prod * (x-dati_x(j))/ (dati_x(i)-dati_x(j))
  endif
enddo
lagrange = prod
return
end function lagrange