Numerical computation of the roots of a polynomial having complex coefficients, based on aberth's method.
this routine approximates the roots of the polynomial
p(x)=a(n+1)x^n+a(n)x^(n-1)+...+a(1), a(j)=cr(j)+i ci(j), i**2=-1,
where a(1) and a(n+1) are nonzero.
the coefficients are complex numbers. the routine is fast, robust against overflow, and allows to deal with polynomials of any degree. overflow situations are very unlikely and may occurr if there exist simultaneously coefficients of moduli close to big and close to small, i.e., the greatest and the smallest positive real(wp) numbers, respectively. in this limit situation the program outputs a warning message. the computation can be speeded up by performing some side computations in single precision, thus slightly reducing the robustness of the program (see the comments in the routine aberth). besides a set of approximations to the roots, the program delivers a set of a-posteriori error bounds which are guaranteed in the most part of cases. in the situation where underflow does not allow to compute a guaranteed bound, the program outputs a warning message and sets the bound to 0. in the situation where the root cannot be represented as a complex(wp) number the error bound is set to -1.
the computation is performed by means of aberth's method according to the formula
x(i)=x(i)-newt/(1-newt*abcorr), i=1,...,n (1)
where newt=p(x(i))/p'(x(i)) is the newton correction and abcorr=
=1/(x(i)-x(1))+...+1/(x(i)-x(i-1))+1/(x(i)-x(i+1))+...+1/(x(i)-x(n))
is the aberth correction to the newton method.
the value of the newton correction is computed by means of the synthetic division algorithm (ruffini-horner's rule) if |x|<=1, otherwise the following more robust (with respect to overflow) formula is applied:
newt=1/(n*y-y**2 r'(y)/r(y)) (2)
where
y=1/x
r(y)=a(1)*y**n+...+a(n)*y+a(n+1) (2')
this computation is performed by the routine newton.
the starting approximations are complex numbers that are
equispaced on circles of suitable radii. the radius of each
circle, as well as the number of roots on each circle and the
number of circles, is determined by applying rouche's theorem
to the functions a(k+1)*x**k and p(x)-a(k+1)*x**k, k=0,...,n.
this computation is performed by the routine start.
if the condition
|p(x(j))|<eps s(|x(j)|) (3)
is satisfied, where s(x)=s(1)+x*s(2)+...+x**n * s(n+1),
s(i)=|a(i)|*(1+3.8*(i-1)), eps is the machine precision (eps=2**-53
for the ieee arithmetic), then the approximation x(j) is not updated
and the subsequent iterations (1) for i=j are skipped.
the program stops if the condition (3) is satisfied for j=1,...,n,
or if the maximum number nitmax of iterations has been reached.
the condition (3) is motivated by a backward rounding error analysis
of the ruffini-horner rule, moreover the condition (3) guarantees
that the computed approximation x(j) is an exact root of a slightly
perturbed polynomial.
for each approximation x of a root, an a-posteriori absolute error
bound r is computed according to the formula
r=n(|p(x)|+eps s(|x|))/|p'(x)| (4)
this provides an inclusion disk of center x and radius r containing a
root.
| Type | Intent | Optional | Attributes | Name | ||
|---|---|---|---|---|---|---|
| integer, | intent(in) | :: | n |
degree of the polynomial. |
||
| complex(kind=wp), | intent(in) | :: | poly(n+1) |
complex vector of n+1 components, |
||
| integer, | intent(in) | :: | nitmax |
the max number of allowed iterations. |
||
| complex(kind=wp), | intent(out) | :: | root(n) |
complex vector of |
||
| real(kind=wp), | intent(out) | :: | radius(n) |
real vector of |
||
| logical, | intent(out) | :: | err(n) |
vector of
the vector |
subroutine polzeros(n, poly, nitmax, root, radius, err) implicit none integer,intent(in) :: n !! degree of the polynomial. complex(wp),intent(in) :: poly(n + 1) !! complex vector of n+1 components, `poly(i)` is the !! coefficient of `x**(i-1), i=1,...,n+1` of the polynomial `p(x)` integer,intent(in) :: nitmax !! the max number of allowed iterations. complex(wp),intent(out) :: root(n) !! complex vector of `n` components, containing the !! approximations to the roots of `p(x)`. real(wp),intent(out) :: radius(n) !! real vector of `n` components, containing the error bounds to !! the approximations of the roots, i.e. the disk of center !! `root(i)` and radius `radius(i)` contains a root of `p(x)`, for !! `i=1,...,n`. `radius(i)` is set to -1 if the corresponding root !! cannot be represented as floating point due to overflow or !! underflow. logical,intent(out) :: err(n) !! vector of `n` components detecting an error condition: !! !! * `err(j)=.true.` if after `nitmax` iterations the stop condition !! (3) is not satisfied for x(j)=root(j); !! * `err(j)=.false.` otherwise, i.e., the root is reliable, !! i.e., it can be viewed as an exact root of a !! slightly perturbed polynomial. !! !! the vector `err` is used also in the routine convex hull for !! storing the abscissae of the vertices of the convex hull. integer :: iter !! number of iterations peformed real(wp) :: apoly(n + 1) !! auxiliary variable: real vector of n+1 components used to store the moduli of !! the coefficients of p(x) and the coefficients of s(x) used !! to test the stop condition (3). real(wp) :: apolyr(n + 1) !! auxiliary variable: real vector of n+1 components used to test the stop !! condition integer :: i, nzeros complex(wp) :: corr, abcorr real(wp) :: amax real(wp),parameter :: eps = epsilon(1.0_wp) real(wp),parameter :: small = tiny(1.0_wp) real(wp),parameter :: big = huge(1.0_wp) ! check consistency of data if (abs(poly(n + 1)) == 0.0_wp) then error stop 'inconsistent data: the leading coefficient is zero' end if if (abs(poly(1)) == 0.0_wp) then error stop 'the constant term is zero: deflate the polynomial' end if ! compute the moduli of the coefficients amax = 0.0_wp do i = 1, n + 1 apoly(i) = abs(poly(i)) amax = max(amax, apoly(i)) apolyr(i) = apoly(i) end do if ((amax) >= (big/(n + 1))) then write (*, *) 'warning: coefficients too big, overflow is likely' end if ! initialize do i = 1, n radius(i) = 0.0_wp err(i) = .true. end do ! select the starting points call start(n, apolyr, root, radius, nzeros, small, big) ! compute the coefficients of the backward-error polynomial do i = 1, n + 1 apolyr(n - i + 2) = eps*apoly(i)*(3.8_wp*(n - i + 1) + 1) apoly(i) = eps*apoly(i)*(3.8_wp*(i - 1) + 1) end do if ((apoly(1) == 0.0_wp) .or. (apoly(n + 1) == 0.0_wp)) then write (*, *) 'warning: the computation of some inclusion radius' write (*, *) 'may fail. this is reported by radius=0' end if do i = 1, n err(i) = .true. if (radius(i) == -1) err(i) = .false. end do ! starts aberth's iterations do iter = 1, nitmax do i = 1, n if (err(i)) then call newton(n, poly, apoly, apolyr, root(i), small, radius(i), corr, err(i)) if (err(i)) then call aberth(n, i, root, abcorr) root(i) = root(i) - corr/(1 - corr*abcorr) else nzeros = nzeros + 1 if (nzeros == n) return end if end if end do end do end subroutine polzeros