program zeri
  ! The program computes the zeros of a polynomial
  implicit none

  real(16) :: x,a,b,x0
  real(16), dimension(:), allocatable :: p
  integer  :: n,k,ios,selettore

  character(7) :: nomefile="pol.dat"

  ! Apro il file pol.dat da cui leggo i dati
  open(unit=8,iostat=ios,file=nomefile,action="read")
  if (ios/=0) then
     write(*,'(a)') "Errore. Impossibile aprire il file ",nomefile
     stop
  end if

  ! Leggo il grado del polinomio e i suoi n+1 coefficienti
  read (8,*) n

  allocate (p(n+1))
  do k=1,n+1
     read (8,*) p(k)

  end do

  ! selettore sceglie l'esame da eseguire
  write(*,'(a)') "Scegli il metodo"
  write(*,'(a)') "1. Metodo della bisezione"
  write(*,'(a)') "2. Metodo delle secanti"
  write(*,'(a)') "3. Metodo di Newton"
  write(*,'(a)') "4. Metodo di Halley"
  read (*,*) selettore

  select case (selettore)

  case(1)
     write(*,'(a)',advance='no') "Inserisci gli estremi dell'intervallo "
     read (*,*) a,b
     call bisezione (n,p,a,b,x)

  case(2)
     write(*,'(a)',advance='no') "Inserisci gli estremi dell'intervallo "
     read (*,*) a,b
     call secanti (n,p,a,b,x)

  case(3)
     write(*,'(a)',advance='no') "Inserisci l'approssimazione iniziale "
     read (*,*) x0
     call newton (n,p,x0)
 
  case(4)
     write(*,'(a)',advance='no') "Inserisci l'approssimazione iniziale "
     read (*,*) x0
     call halley (n,p,x0)

  end select
end program zeri

subroutine bisezione (n,p,a,b,x)
  implicit none

  integer  :: n,iter
  real(16) :: a,b,c,x,eps=1.q-32,fa,fb,fc,horner
  real(16),dimension(n) :: p

  fa=horner(n,p,a);fb=horner(n,p,b)

  if (fa*fb>0) then
     write(*,'(a,f7.3,a,f7.3,a)') "La funzione non ha zeri nell'intervallo [",a,",",b,"]"
     stop
  end if

  iter=0
  do
     if (abs(b-a)<eps) exit
     iter=iter+1
     c=(a+b)/2
     fc=horner(n,p,c)
     if (fa*fc>=0) then
        a=c;fa=fc
     else
        b=c;fb=fc
     end if
!     ! Stampo su schermo l'approssimazione al passo iter
     write (*,'(a,i3,a,f)') "Iterazione = ",iter," Radice = ",c
  end do

  x=c
  write (*,*) "Risultato con il metodo di bisezione"
  write (*,*) "Iterazioni = ",iter," Radice = ",x
end subroutine bisezione

subroutine secanti (n,p,a,b,x)
  implicit none

  integer  :: n,iter
  real(16) :: a,b,c,x,eps=1.q-32,fa,fb,fc,horner
  real(16),dimension(n) :: p

  fa=horner(n,p,a);fb=horner(n,p,b)

  if (fa*fb>0) then
     write(*,'(a,f7.3,a,f7.3,a)') "La funzione non ha zeri nell'intervallo [",a,",",b,"]"
     stop
  end if

  iter=0
  do
     if (abs(b-a)<eps) exit
     iter=iter+1
     c=b-(b-a)*fb/(fb-fa)
     fc=horner(n,p,c)
     a=b;fa=fb
     b=c;fb=fc
!     ! Stampo su schermo l'approssimazione al passo iter
     write (*,'(a,i3,a,f)') "Iterazione = ",iter," Radice = ",c
  end do

  x=c
  write (*,*) "Risultato con il metodo delle secanti"
  write (*,*) "Iterazioni = ",iter," Radice = ",x
end subroutine secanti

subroutine newton (n,p,x0)
  implicit none
  
  integer :: n,iter,k
  real(16),dimension (n+1) :: p,b,c
  real(16) :: x0,xold,xnew,eps=1.q-32

  iter=0;xnew=x0;xold=x0+1
  do
     if (abs(xold-xnew)<eps) exit
     xold=xnew
     iter=iter+1

     ! Calcolo f(xold) e f'(xold) con un simil-horner
     b(n+1)=p(n+1);c(n+1)=b(n+1)
     horner: do k=n,2,-1
        b(k)=b(k+1)*xold+p(k)
!       c(k)=b(k+1)*xold+b(k) Funzione errata che converge :) Spiegate perche'
        c(k)=c(k+1)*xold+b(k)
     end do horner
     b(1)=b(2)*xold+p(1)

     xnew=xold-b(1)/c(2)
!     ! Stampo su schermo l'approssimazione al passo iter
     write (*,'(a,i3,a,f)') "Iterazione = ",iter," Radice = ",xnew
  end do

end subroutine newton

subroutine halley (n,p,x0)
  implicit none
  
  integer :: n,iter,k
  real(16),dimension (n+1) :: p,b,c,d
  real(16) :: x0,xold,xnew,eps=1.q-32

  iter=0;xnew=x0;xold=x0+1
  do
     if (abs(xold-xnew)<eps) exit
     xold=xnew
     iter=iter+1

     ! Calcolo f(xold) e f'(xold) con un simil-horner
     b(n+1)=p(n+1);c(n+1)=b(n+1)
     horner: do k=n,3,-1
        b(k)=b(k+1)*xold+p(k)
        c(k)=c(k+1)*xold+b(k)
        d(k)=d(k+1)*xold+c(k)
     end do horner
     b(2)=b(3)*xold+p(2)
     c(2)=c(3)*xold+b(2)
     b(1)=b(2)*xold+p(1)

    xnew=xold-(b(1)/c(2))/(1-b(1)*c(3)/c(2)**2)
    ! Stampo su schermo l'approssimazione al passo iter
     write (*,'(a,i3,a,f)') "Iterazione = ",iter," Radice = ",xnew
  end do

end subroutine halley

function horner(n,p,x)
  ! Calcola il valore p(x) mediante il metodo di Horner
  implicit none
  
  integer  :: n,k
  real(16) :: x,horner
  real(16),dimension(n+1):: p

  horner=p(n+1)
  do k=n,1,-1
     horner=horner*x+p(k)
  end do

end function horner

