newton Subroutine

    subroutine newton(n, poly, apoly, apolyr, z, small, radius, corr, again)

        !! compute the newton's correction, the inclusion radius (4) and checks
        !! the stop condition (3)

        implicit none

        integer,intent(in) :: n !! degree of the polynomial p(x)
        complex(wp),intent(in) :: poly(n + 1) !! coefficients of the polynomial p(x)
        real(wp),intent(in) :: apoly(n + 1) !! upper bounds on the backward perturbations on the
                                            !! coefficients of p(x) when applying ruffini-horner's rule
        real(wp),intent(in) :: apolyr(n + 1) !! upper bounds on the backward perturbations on the
                                             !! coefficients of p(x) when applying (2), (2')
        complex(wp),intent(in) :: z !! value at which the newton correction is computed
        real(wp),intent(in) :: small !! the min positive real(wp), small=2**(-1074) for the ieee.
        real(wp),intent(out) :: radius !! upper bound to the distance of z from the closest root of
                                       !! the polynomial computed according to (4).
        complex(wp),intent(out) :: corr !! newton's correction
        logical,intent(out) :: again !! this variable is .true. if the computed value p(z) is
                                     !! reliable, i.e., (3) is not satisfied in z. again is
                                     !! .false., otherwise.

        integer :: i
        complex(wp) :: p, p1, zi, den, ppsp
        real(wp) :: ap, az, azi, absp

        az = abs(z)
        ! if |z|<=1 then apply ruffini-horner's rule for p(z)/p'(z)
        ! and for the computation of the inclusion radius
        if (az <= 1.0_wp) then
            p = poly(n + 1)
            ap = apoly(n + 1)
            p1 = p
            do i = n, 2, -1
                p = p*z + poly(i)
                p1 = p1*z + p
                ap = ap*az + apoly(i)
            end do
            p = p*z + poly(1)
            ap = ap*az + apoly(1)
            corr = p/p1
            absp = abs(p)
            ap = ap
            again = (absp > (small + ap))
            if (.not. again) radius = n*(absp + ap)/abs(p1)
        else
            ! if |z|>1 then apply ruffini-horner's rule to the reversed polynomial
            ! and use formula (2) for p(z)/p'(z). analogously do for the inclusion
            ! radius.
            zi = 1.0_wp/z
            azi = 1.0_wp/az
            p = poly(1)
            p1 = p
            ap = apolyr(n + 1)
            do i = n, 2, -1
                p = p*zi + poly(n - i + 2)
                p1 = p1*zi + p
                ap = ap*azi + apolyr(i)
            end do
            p = p*zi + poly(n + 1)
            ap = ap*azi + apolyr(1)
            absp = abs(p)
            again = (absp > (small + ap))
            ppsp = (p*z)/p1
            den = n*ppsp - 1
            corr = z*(ppsp/den)
            if (again) return
            radius = abs(ppsp) + (ap*az)/abs(p1)
            radius = n*radius/abs(den)
            radius = radius*az
        end if

    end subroutine newton