aberth Subroutine

    subroutine aberth(n, j, root, abcorr)

        !! compute the aberth correction. to save time, the reciprocation of
        !! root(j)-root(i) could be performed in single precision (complex*8)
        !! in principle this might cause overflow if both root(j) and root(i)
        !! have too small moduli.

        implicit none

        integer,intent(in) :: n !! degree of the polynomial
        integer,intent(in) :: j  !! index of the component of root with respect to which the
                                 !! aberth correction is computed
        complex(wp),intent(in) :: root(n) !! vector containing the current approximations to the roots
        complex(wp),intent(out) :: abcorr !! aberth's correction (compare (1))

        integer :: i
        complex(wp) :: z, zj

        abcorr = 0.0_wp
        zj = root(j)
        do i = 1, j - 1
            z = zj - root(i)
            abcorr = abcorr + 1.0_wp/z
        end do
        do i = j + 1, n
            z = zj - root(i)
            abcorr = abcorr + 1.0_wp/z
        end do

    end subroutine aberth