MODULE poly_zeroes USE quadruple_precision USE quad_prec_complex IMPLICIT NONE PRIVATE :: aberth, newton, start, cnvex, cmerge, left, right, ctest PUBLIC :: polzeros CONTAINS !************************************************************************ ! NUMERICAL COMPUTATION OF THE ROOTS OF A POLYNOMIAL HAVING * ! COMPLEX COEFFICIENTS, BASED ON ABERTH'S METHOD. * ! Version 1.4, June 1996 * ! (D. Bini, Dipartimento di Matematica, Universita' di Pisa) * ! (bini@dm.unipi.it) * ! ! *************************************************************************** ! * All the software contained in this library is protected by copyright. * ! * Permission to use, copy, modify, and distribute this software for any * ! * purpose without fee is hereby granted, provided that this entire notice * ! * is included in all copies of any software which is or includes a copy * ! * or modification of this software and in all copies of the supporting * ! * documentation for such software. * ! *************************************************************************** ! * THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR IMPLIED * ! * WARRANTY. IN NO EVENT, NEITHER THE AUTHORS, NOR THE PUBLISHER, NOR ANY * ! * MEMBER OF THE EDITORIAL BOARD OF THE JOURNAL "NUMERICAL ALGORITHMS", * ! * NOR ITS EDITOR-IN-CHIEF, BE LIABLE FOR ANY ERROR IN THE SOFTWARE, ANY * ! * MISUSE OF IT OR ANY DAMAGE ARISING OUT OF ITS USE. THE ENTIRE RISK OF * ! * USING THE SOFTWARE LIES WITH THE PARTY DOING SO. * ! *************************************************************************** ! * ANY USE OF THE SOFTWARE CONSTITUTES ACCEPTANCE OF THE TERMS OF THE * ! * ABOVE STATEMENT. * ! *************************************************************************** ! ! AUTHOR: ! ! DARIO ANDREA BINI ! UNIVERSITY OF PISA, ITALY ! E-MAIL: bini@dm.unipi.it ! ! REFERENCE: ! ! - NUMERICAL COMPUTATION OF POLYNOMIAL ZEROS BY MEANS OF ABERTH'S METHOD ! NUMERICAL ALGORITHMS, 13 (1996), PP. 179-200 ! ! Work performed under the support of the ESPRIT BRA project 6846 POSSO * ! ************************************************************************* ! ! This version, which is compatible with Lahey's free ELF90 compiler, ! is by Alan Miller, CSIRO Mathematical & Information Sciences, ! Private Bag 10, Clayton South MDC, Victoria, Australia 3169. ! Alan.Miller @ mel.dms.csiro.au http://www.ozemail.com.au/~milleraj ! Latest revision of ELF90 version - 19 January 1998 ! ! ************************************************************************* ! ! QUADRUPLE - PRECISION VERSION ! !************************************************************************ !********** SUBROUTINES AND FUNCTIONS *********** !************************************************************************ ! The following routines/functions are listed: * ! POLZEROS : computes polynomial roots by means of Aberth's method * ! ABERTH : computes the Aberth correction * ! NEWTON : computes p(x)/p'(x) by means of Ruffini-Horner's rule * ! START : Selects N starting points by means of Rouche's theorem * ! CNVEX : Computes the convex hull, used by START * ! CMERGE : Used by CNVEX * ! LEFT : Used by CMERGE * ! RIGHT : Used by CMERGE * ! CTEST : Convexity test, Used by CMERGE * !************************************************************************ ! * ! * !************************************************************************ !********************** SUBROUTINE POLZEROS ***************************** !************************************************************************ ! GENERAL COMMENTS * !************************************************************************ ! 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*16 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*8 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*16 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. * !************************************************************************ ! STOP CONDITION * !************************************************************************ ! If the condition * ! |p(x(j))| amax) amax = apoly(i) apolyr(i) = apoly(i) END DO IF (amax%hi >= big/(n+1)) THEN WRITE(*,*)'WARNING: COEFFICIENTS TOO BIG, OVERFLOW IS LIKELY' WRITE(2,*)'WARNING: COEFFICIENTS TOO BIG, OVERFLOW IS LIKELY' END IF ! Initialize DO i = 1, n radius(i) = zero err(i) = .true. END DO ! Select the starting points CALL start(n, apolyr, root, radius, nzeros, small, big, err) ! Compute the coefficients of the backward-error polynomial DO i = 1, n+1 apolyr(n-i+2) = eps*apoly(i)*(3.8_dp*(n-i+1) + 1._dp) apoly(i) = eps*apoly(i)*(3.8_dp*(i-1) + 1._dp) END DO IF ((apoly(1)%hi == 0.0_dp) .OR. (apoly(n+1)%hi == 0.0_dp)) THEN WRITE(*,*)'WARNING: THE COMPUTATION OF SOME INCLUSION RADIUS' WRITE(*,*)'MAY FAIL. THIS IS REPORTED BY RADIUS = 0' WRITE(2,*)'WARNING: THE COMPUTATION OF SOME INCLUSION RADIUS' WRITE(2,*)'MAY FAIL. THIS IS REPORTED BY RADIUS = 0' END IF DO i = 1,n err(i) = .true. IF (ABS( radius(i)%hi + 1.0_dp) < eps) 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 / (qc(one,zero) - corr*abcorr) ELSE nzeros = nzeros + 1 IF (nzeros == n) RETURN END IF END IF END DO END DO RETURN END SUBROUTINE polzeros !************************************************************************ ! SUBROUTINE NEWTON * !************************************************************************ ! Compute the Newton's correction, the inclusion radius (4) and checks * ! the stop condition (3) * !************************************************************************ ! Input variables: * ! N : degree of the polynomial p(x) * ! POLY : coefficients of the polynomial p(x) * ! APOLY : upper bounds on the backward perturbations on the * ! coefficients of p(x) when applying Ruffini-Horner's rule * ! APOLYR: upper bounds on the backward perturbations on the * ! coefficients of p(x) when applying (2), (2') * ! Z : value at which the Newton correction is computed * ! SMALL : the min positive TYPE (quad) ::, SMALL=2**(-1074) for the IEEE. * !************************************************************************ ! Output variables: * ! RADIUS: upper bound to the distance of Z from the closest root of * ! the polynomial computed according to (4). * ! CORR : Newton's correction * ! 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. * !************************************************************************ SUBROUTINE newton(n, poly, apoly, apolyr, z, small, radius, corr, again) INTEGER, INTENT(IN) :: n TYPE (qc), INTENT(IN) :: poly(:), z TYPE (qc), INTENT(OUT) :: corr TYPE (quad), INTENT(IN) :: apoly(:), apolyr(:) REAL (dp), INTENT(IN) :: small TYPE (quad), INTENT(OUT) :: radius LOGICAL, INTENT(OUT) :: again ! Local variables INTEGER :: i TYPE (qc) :: p, p1, zi, den, ppsp TYPE (quad) :: ap, az, azi, absp TYPE (quad), PARAMETER :: zero = quad( 0.0_dp, 0.0_dp ), & one = quad( 1.0_dp, 0.0_dp ) 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%hi <= 1.0_dp) 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 = ABS(ap) again = (absp > quad(small, 0.0_dp) + ap) IF (.NOT.again) radius = DBLE(n)*(absp + ap) / ABS(p1) RETURN 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 = qc(one,zero) / z azi = one / 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 > quad(small, 0.0_dp) + ap) ppsp = (p*z)/p1 den = qc(quad(DBLE(n),0.0_dp), zero) * ppsp - qc(one,zero) corr = z*(ppsp/den) IF (again)RETURN radius = ABS(ppsp) + (ap*az) / ABS(p1) radius = n*radius/ABS(den) radius = radius*az END IF RETURN END SUBROUTINE newton !************************************************************************ ! SUBROUTINE ABERTH * !************************************************************************ ! 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. * !************************************************************************ ! Input variables: * ! N : degree of the polynomial * ! ROOT : vector containing the current approximations to the roots * ! J : index of the component of ROOT with respect to which the * ! Aberth correction is computed * !************************************************************************ ! Output variable: * ! ABCORR: Aberth's correction (compare (1)) * !************************************************************************ SUBROUTINE aberth(n, j, root, abcorr) INTEGER, INTENT(IN) :: n, j TYPE (qc), INTENT(IN) :: root(:) TYPE (qc), INTENT(OUT) :: abcorr ! Local variables INTEGER :: i TYPE (qc) :: z, zj TYPE (quad), PARAMETER :: zero = quad( 0.0_dp, 0.0_dp ), & one = quad( 1.0_dp, 0.0_dp ) abcorr = CMPLX(zero, zero) zj = root(j) DO i = 1, j-1 z = zj - root(i) abcorr = abcorr + qc(one,zero) / z END DO DO i = j+1, n z = zj - root(i) abcorr = abcorr + qc(one,zero) / z END DO RETURN END SUBROUTINE aberth !************************************************************************ ! SUBROUTINE START * !************************************************************************ ! Compute the starting approximations of the roots * !************************************************************************ ! Input variables: * ! N : number of the coefficients of the polynomial * ! A : moduli of the coefficients of the polynomial * ! SMALL : the min positive REAL (dp), SMALL=2**(-1074) for the IEEE. * * ! BIG : the max REAL (dp), BIG=2**1023 for the IEEE standard. * ! Output variables: * ! Y : starting approximations * ! RADIUS: if a component is -1 then the corresponding root has a * ! too big or too small modulus in order to be represented * ! as double float with no overflow/underflow * ! NZ : number of roots which cannot be represented without * ! overflow/underflow * ! Auxiliary variables: * ! H : needed for the computation of the convex hull * !************************************************************************ ! This routines selects starting approximations along circles center at * ! 0 and having suitable radii. The computation of the number of circles * ! and of the corresponding radii is performed by computing the upper * ! convex hull of the set (i,log(A(i))), i=1,...,n+1. * !************************************************************************ SUBROUTINE start(n, a, y, radius, nz, small, big, h) INTEGER, INTENT(IN) :: n INTEGER, INTENT(OUT) :: nz LOGICAL, INTENT(OUT) :: h(:) TYPE (qc), INTENT(OUT) :: y(:) REAL (dp), INTENT(IN) :: small, big TYPE (quad), INTENT(IN OUT) :: a(:) TYPE (quad), INTENT(OUT) :: radius(:) ! Local variables INTEGER :: i, iold, nzeros, j, jj TYPE (quad) :: r, th, ang, temp REAL (dp) :: xsmall, xbig TYPE (quad), PARAMETER :: pi2 = twopi, sigma = quad( 0.7_dp, 0.0_dp ) xsmall = LOG(small) xbig = LOG(big) nz = 0 ! Compute the logarithm A(I) of the moduli of the coefficients of ! the polynomial and then the upper convex hull of the set (A(I),I) DO i = 1, n+1 IF (a(i)%hi /= 0.0_dp) THEN a(i) = LOG(a(i)) ELSE a(i) = quad(-1.e50_dp, 0.0_dp) END IF END DO CALL cnvex(n+1, a, h) ! Given the upper convex hull of the set (A(I),I) compute the moduli ! of the starting approximations by means of Rouche's theorem iold = 1 th = pi2 / DBLE(n) DO i = 2, n+1 IF (h(i)) THEN nzeros = i - iold temp = (a(iold) - a(i)) / DBLE(nzeros) ! Check if the modulus is too small IF ((temp%hi < -xbig) .AND. (temp%hi >= xsmall)) THEN WRITE(*,*)'WARNING:',nzeros,' ZERO(S) ARE TOO SMALL TO' WRITE(*,*)'REPRESENT THEIR INVERSES AS TYPE (qc), THEY' WRITE(*,*)'ARE REPLACED BY SMALL NUMBERS, THE CORRESPONDING' WRITE(*,*)'RADII ARE SET TO -1' WRITE(2,*)'WARNING:',nzeros,' ZERO(S) ARE TOO SMALL TO ' WRITE(2,*)'REPRESENT THEIR INVERSES AS TYPE (qc), THEY' WRITE(2,*)'ARE REPLACED BY SMALL NUMBERS, THE CORRESPONDING' WRITE(2,*)'RADII ARE SET TO -1' nz = nz + nzeros r = quad(1.0_dp, 0.0_dp) / big END IF IF (temp%hi < xsmall) THEN nz = nz + nzeros WRITE(*,*)'WARNING: ',nzeros,' ZERO(S) ARE TOO SMALL TO BE' WRITE(*,*)'REPRESENTED AS TYPE (qc), THEY ARE SET TO 0' WRITE(*,*)'THE CORRESPONDING RADII ARE SET TO -1' WRITE(2,*)'WARNING: ',nzeros,' ZERO(S) ARE TOO SMALL TO BE' WRITE(2,*)'REPRESENTED AS TYPE (qc), THEY ARE SET 0' WRITE(2,*)'THE CORRESPONDING RADII ARE SET TO -1' END IF ! Check if the modulus is too big IF (temp%hi > xbig) THEN r = quad(big, 0.0_dp) nz = nz + nzeros WRITE(*,*)'WARNING: ', nzeros, ' ZEROS(S) ARE TOO BIG TO BE' WRITE(*,*)'REPRESENTED AS TYPE (qc),' WRITE(*,*)'THE CORRESPONDING RADII ARE SET TO -1' WRITE(2,*)'WARNING: ',nzeros, ' ZERO(S) ARE TOO BIG TO BE' WRITE(2,*)'REPRESENTED AS TYPE (qc),' WRITE(2,*)'THE CORRESPONDING RADII ARE SET TO -1' END IF IF ((temp%hi <= xbig) .AND. (temp%hi > MAX(-xbig, xsmall))) THEN r = EXP(temp) END IF ! Compute NZEROS approximations equally distributed in the disk of radius R ang = pi2 / DBLE(nzeros) DO j = iold, i-1 jj = j - iold + 1 IF (r%hi <= 1.0_dp/big .OR. r == quad(big, 0.0_dp)) & radius(j) = quad(-1.0_dp, 0.0_dp) y(j)%qr = r * COS(ang*jj + th*i + sigma) y(j)%qi = r * SIN(ang*jj + th*i + sigma) END DO iold = i END IF END DO RETURN END SUBROUTINE start !************************************************************************ ! SUBROUTINE CNVEX * !************************************************************************ ! Compute the upper convex hull of the set (i,a(i)), i.e., the set of * ! vertices (i_k,a(i_k)), k=1,2,...,m, such that the points (i,a(i)) lie * ! below the straight lines passing through two consecutive vertices. * ! The abscissae of the vertices of the convex hull equal the indices of * ! the TRUE components of the logical output vector H. * ! The used method requires O(nlog n) comparisons and is based on a * ! divide-and-conquer technique. Once the upper convex hull of two * ! contiguous sets (say, {(1,a(1)),(2,a(2)),...,(k,a(k))} and * ! {(k,a(k)), (k+1,a(k+1)),...,(q,a(q))}) have been computed, then * ! the upper convex hull of their union is provided by the subroutine * ! CMERGE. The program starts with sets made up by two consecutive * ! points, which trivially constitute a convex hull, then obtains sets * ! of 3,5,9... points, up to arrive at the entire set. * ! The program uses the subroutine CMERGE; the subroutine CMERGE uses * ! the subroutines LEFT, RIGHT and CTEST. The latter tests the convexity * ! of the angle formed by the points (i,a(i)), (j,a(j)), (k,a(k)) in the * ! vertex (j,a(j)) up to within a given tolerance TOLER, where iI and H(IL) is TRUE. * !************************************************************************ !************************************************************************ ! Input variables: * ! N : length of the vector H * ! H : vector of logical * ! I : integer * !************************************************************************ ! Output variable: * ! IR : minimum integer such that IR>I, H(IR)=.TRUE. * !************************************************************************ SUBROUTINE right(n, h, i, ir) INTEGER, INTENT(IN) :: n, i INTEGER, INTENT(OUT) :: ir LOGICAL, INTENT(IN) :: h(:) DO ir = i+1, n IF (h(ir)) RETURN END DO RETURN END SUBROUTINE right !************************************************************************ ! SUBROUTINE CMERGE * !************************************************************************ ! Given the upper convex hulls of two consecutive sets of pairs * ! (j,A(j)), compute the upper convex hull of their union * !************************************************************************ ! Input variables: * ! N : length of the vector A * ! A : vector defining the points (j,A(j)) * ! I : abscissa of the common vertex of the two sets * ! M : the number of elements of each set is M+1 * !************************************************************************ ! Input/Output variable: * ! H : vector defining the vertices of the convex hull, i.e., * ! H(j) is .TRUE. if (j,A(j)) is a vertex of the convex hull * ! This vector is used also as output. * !************************************************************************ SUBROUTINE cmerge(n, a, i, m, h) INTEGER, INTENT(IN) :: n, m, i LOGICAL, INTENT(IN OUT) :: h(:) TYPE (quad), INTENT(IN) :: a(:) ! Local variables INTEGER :: ir, il, irr, ill LOGICAL :: tstl, tstr ! at the left and the right of the common vertex (I,A(I)) determine ! the abscissae IL,IR, of the closest vertices of the upper convex ! hull of the left and right sets, respectively CALL left(h, i, il) CALL right(n, h, i, ir) ! check the convexity of the angle formed by IL,I,IR IF (ctest(a, il, i, ir)) THEN RETURN ELSE ! continue the search of a pair of vertices in the left and right ! sets which yield the upper convex hull h(i) = .false. DO IF (il == (i-m)) THEN tstl = .true. ELSE CALL left(h, il, ill) tstl = ctest(a, ill, il, ir) END IF IF (ir == MIN(n, i+m)) THEN tstr = .true. ELSE CALL right(n, h, ir, irr) tstr = ctest(a, il, ir, irr) END IF h(il) = tstl h(ir) = tstr IF (tstl.AND.tstr) RETURN IF (.NOT.tstl) il = ill IF (.NOT.tstr) ir = irr END DO END IF RETURN END SUBROUTINE cmerge !************************************************************************ ! FUNCTION CTEST * !************************************************************************ ! Test the convexity of the angle formed by (IL,A(IL)), (I,A(I)), * ! (IR,A(IR)) at the vertex (I,A(I)), up to within the tolerance * ! TOLER. If convexity holds then the function is set to .TRUE., * ! otherwise CTEST=.FALSE. The parameter TOLER is set to 0.4 by default. * !************************************************************************ ! Input variables: * ! A : vector of double * ! IL,I,IR : integers such that IL < I < IR * !************************************************************************ ! Output: * ! .TRUE. if the angle formed by (IL,A(IL)), (I,A(I)), (IR,A(IR)) at * ! the vertex (I,A(I)), is convex up to within the tolerance * ! TOLER, i.e., if * ! (A(I)-A(IL))*(IR-I) - (A(IR)-A(I))*(I-IL) > TOLER. * ! .FALSE., otherwise. * !************************************************************************ FUNCTION ctest(a, il, i, ir) RESULT(OK) INTEGER, INTENT(IN) :: i, il, ir TYPE (quad), INTENT(IN) :: a(:) LOGICAL :: OK ! Local variables TYPE (quad) :: s1, s2 REAL (dp), PARAMETER :: toler = 0.4_dp s1 = a(i) - a(il) s2 = a(ir) - a(i) s1 = s1*(ir-i) s2 = s2*(i-il) OK = .false. IF (s1 > (s2 + toler)) OK = .true. RETURN END FUNCTION ctest END MODULE poly_zeroes