program rootfinding
  ! Calcolo simultaneo delle radici di un polinomio
  implicit none

  integer :: n,i
  complex(8),allocatable,dimension(:) :: p,x

  ! In questa versione viene richiesto il grado del polinomio
  ! Il polinomio e` x^n-1
!  write(*,'(a)') "Inerisci il grado del polinomio "
!  read(*,*) n
  n=8
  allocate(p(n+1),x(n))
  p=0;p(1)=-1;p(n+1)=1

  call aberth(n,p,x)
  do i=1,n
     write(*,*) x(i)
  end do
end program rootfinding

subroutine aberth (n,p,xnew)
  ! Calcola le radici x di un polinomio complesso p di grado n
  ! tramite il metodo di Aberth
  implicit none

  integer :: n,k,iter
  complex(8),dimension(n+1) :: p
  complex(8),dimension(n)   :: xnew,xold
  complex(8),parameter :: i=(0.d0,1.d0)
  complex(8) :: correzione
  real(8),parameter :: pi=3.141592653589793d0

  do k=1,n
     xold(k)=exp((2*i*pi*k+0.1)/n)
  end do
  ! xold(1:n)=exp((2*i*pi*(1:n)+0.1)/n)
  iter=0
  do 
     do k=1,n
        xnew(k)=xold(k)-correzione(n,k,p,xold)
     end do
     iter=iter+1
     if (maxval(abs(xnew-xold))<1.0d-12) exit
     if (iter >= 100) exit
     xold=xnew
  end do
end subroutine aberth

function correzione (n,i,p,xk)
  ! Calcola la correzione di Aberth
  implicit none

  integer :: n,i
  complex(8) :: x,y,z,frac,sum,correzione
  complex(8),dimension(n)   :: xk
  complex(8),dimension(n+1) :: p

  integer :: h,j

  x=xk(i)
  ! Calcolo y=p(x)
  y=p(n+1)
  do h=n,1,-1
     y=y*x+p(h)
  end do
  ! Calcolo z=p'(x)
  z=p(n+1)*n
  do h=n,2,-1
     z=z*x+p(h)*(h-1)
  end do
  ! Calcolo frac=p'/p
  frac=z/y
  ! Calcolo la sommatoria
  sum=0
  do j=1,n
     if (j/=i) sum=sum+1/(xk(i)-xk(j))
  end do

  correzione=1/(frac-sum)
end function correzione
