MODULE PolynomialRoots
! ---------------------------------------------------------------------------
! PURPOSE - Solve for the roots of a polynomial equation with real
! coefficients, up to quartic order. Returns a code indicating the nature
! of the roots found.
! AUTHORS - Alfred H. Morris, Naval Surface Weapons Center, Dahlgren,VA
! William L. Davis, Naval Surface Weapons Center, Dahlgren,VA
! Alan Miller, CSIRO Mathematical & Information Sciences
! CLAYTON, VICTORIA, AUSTRALIA 3169
! http://www.mel.dms.csiro.au/~alan
! Ralph L. Carmichael, Public Domain Aeronautical Software
! http://www.pdas.com
! REVISION HISTORY
! DATE VERS PERSON STATEMENT OF CHANGES
! ?? 1.0 AHM&WLD Original coding
! 27Feb97 1.1 AM Converted to be compatible with ELF90
! 12Jul98 1.2 RLC Module format; numerous style changes
! 4Jan99 1.3 RLC Made the tests for zero constant term exactly zero
IMPLICIT NONE
INTEGER,PARAMETER,PRIVATE:: SP=selected_real_kind(6), DP=selected_real_kind(12)
REAL(DP),PARAMETER,PRIVATE:: ZERO=0.0D0, FOURTH=0.25D0, HALF=0.5D0
REAL(DP),PARAMETER,PRIVATE:: ONE=1.0D0, TWO=2.0D0, THREE=3.0D0, FOUR=4.0D0
COMPLEX(DP),PARAMETER,PRIVATE:: CZERO=(0.D0,0.D0)
REAL(DP),PARAMETER,PRIVATE:: EPS=EPSILON(ONE)
CHARACTER(LEN=*),PARAMETER,PUBLIC:: POLYROOTS_VERSION= "1.3 (4 Jan 1999)"
INTEGER,PRIVATE:: outputCode
! =0 degenerate equation
! =1 one real root
! =21 two identical real roots
! =22 two distinct real roots
! =23 two complex roots
! =31 multiple real roots
! =32 one real and two complex roots
! =33 three distinct real roots
! =41
! =42 two real and two complex roots
! =43
! =44 four complex roots
PRIVATE:: CubeRoot
PUBLIC:: LinearRoot
PRIVATE:: OneLargeTwoSmall
PUBLIC:: QuadraticRoots
PUBLIC:: CubicRoots
PUBLIC:: QuarticRoots
PUBLIC:: SolvePolynomial
!----------------------------------------------------------------------------
INTERFACE Swap
MODULE PROCEDURE SwapDouble, SwapSingle
END INTERFACE
CONTAINS
!+
FUNCTION CubeRoot(x) RESULT(f)
! ---------------------------------------------------------------------------
! PURPOSE - Compute the Cube Root of a REAL(DP) number. If the argument is
! negative, then the cube root is also negative.
REAL(DP),INTENT(IN) :: x
REAL(DP):: f
!----------------------------------------------------------------------------
IF (x < ZERO) THEN
f=-EXP(LOG(-x)/THREE)
ELSE IF (x > ZERO) THEN
f=EXP(LOG(x)/THREE)
ELSE
f=ZERO
END IF
RETURN
END Function CubeRoot ! ---------------------------------------------------
!+
SUBROUTINE LinearRoot(a, z)
! ---------------------------------------------------------------------------
! PURPOSE - COMPUTES THE ROOTS OF THE REAL POLYNOMIAL
! A(1) + A(2)*Z
! AND STORES THE RESULTS IN Z. It is assumed that a(2) is non-zero.
REAL,INTENT(IN),DIMENSION(:):: a
REAL,INTENT(OUT):: z
!----------------------------------------------------------------------------
IF (a(2)==0.0) THEN
z=0.
ELSE
z=-a(1)/a(2)
END IF
RETURN
END Subroutine LinearRoot ! -----------------------------------------------
!+
SUBROUTINE OneLargeTwoSmall(a1,a2,a4,w, z)
! ---------------------------------------------------------------------------
! PURPOSE - Compute the roots of a cubic when one root, w, is known to be
! much larger in magnitude than the other two
REAL,INTENT(IN):: a1,a2,a4
REAL,INTENT(IN):: w
COMPLEX,INTENT(OUT),DIMENSION(:):: z
REAL(DP),DIMENSION(3):: aq
!----------------------------------------------------------------------------
aq(1)=a1
aq(2)=a2+a1/w
aq(3)=-a4*w
CALL QuadraticRoots(real(aq), z)
z(3)=CMPLX(w,0.)
IF (AIMAG(z(1)) == 0.) RETURN
z(3)=z(2)
z(2)=z(1)
z(1)=CMPLX(w,0.)
RETURN
END Subroutine OneLargeTwoSmall ! -----------------------------------------
!+
SUBROUTINE QuadraticRoots(a, z)
! ---------------------------------------------------------------------------
! PURPOSE - COMPUTES THE ROOTS OF THE REAL POLYNOMIAL
! A(1) + A(2)*Z + A(3)*Z**2
! AND STORES THE RESULTS IN Z. IT IS ASSUMED THAT A(3) IS NONZERO.
REAL,INTENT(IN),DIMENSION(:):: a
COMPLEX,INTENT(OUT),DIMENSION(:):: z
REAL(DP):: d, r, w, x, y
!----------------------------------------------------------------------------
IF(a(1)==0.0) THEN ! EPS is a global module constant (private)
z(1) = CZERO ! one root is obviously zero
z(2) = CMPLX(-a(2)/a(3), 0.) ! remainder is a linear eq.
outputCode=21 ! two identical real roots
RETURN
END IF
d = a(2)*a(2) - FOUR*a(1)*a(3) ! the discriminant
IF (ABS(d) <= TWO*eps*a(2)*a(2)) THEN
z(1) = CMPLX(-0.5*a(2)/a(3), 0.) ! discriminant is tiny
z(2) = z(1)
outputCode=22 ! two distinct real roots
RETURN
END IF
r = SQRT(ABS(d))
IF (d < ZERO) THEN
x = -HALF*a(2)/a(3) ! negative discriminant => roots are complex
y = ABS(HALF*r/a(3))
z(1) = CMPLX(x, y, DP)
z(2) = CMPLX(x,-y, DP) ! its conjugate
outputCode=23 ! COMPLEX ROOTS
RETURN
END IF
IF (a(2) /= ZERO) THEN ! see Numerical Recipes, sec. 5.5
w = -(a(2) + SIGN(r,dble(a(2))))
z(1) = CMPLX(TWO*a(1)/w, ZERO, DP)
z(2) = CMPLX(HALF*w/a(3), ZERO, DP)
outputCode=22 ! two real roots
RETURN
END IF
x = ABS(HALF*r/a(3)) ! a(2)=0 if you get here
z(1) = CMPLX( x, ZERO, DP)
z(2) = CMPLX(-x, ZERO, DP)
outputCode=22
RETURN
END Subroutine QuadraticRoots ! -------------------------------------------
!+
SUBROUTINE CubicRoots(a, z)
!----------------------------------------------------------------------------
! PURPOSE - Compute the roots of the real polynomial
! A(1) + A(2)*Z + A(3)*Z**2 + A(4)*Z**3
REAL,INTENT(IN),DIMENSION(:):: a
COMPLEX,INTENT(OUT),DIMENSION(:):: z
REAL(DP),PARAMETER:: RT3=1.7320508075689D0 ! (Sqrt(3)
REAL:: aq(3)
REAL(DP):: arg, c, cf, d, p, p1, q, q1
REAL(DP):: r, ra, rb, rq, rt
REAL(DP):: r1, s, sf, sq, sum, t, tol, t1, w
REAL(DP):: w1, w2, x, x1, x2, x3, y, y1, y2, y3
! NOTE - It is assumed that a(4) is non-zero. No test is made here.
!----------------------------------------------------------------------------
IF (a(1)==0.0) THEN
z(1) = CZERO ! one root is obviously zero
CALL QuadraticRoots(a(2:4), z(2:3)) ! remaining 2 roots here
RETURN
END IF
p = a(3)/(THREE*a(4))
q = a(2)/a(4)
r = a(1)/a(4)
tol = FOUR*EPS
c = ZERO
t = a(2) - p*a(3)
IF (ABS(t) > tol*ABS(a(2))) c = t/a(4)
t = TWO*p*p - q
IF (ABS(t) <= tol*ABS(q)) t = ZERO
d = r + p*t
IF (ABS(d) <= tol*ABS(r)) GO TO 110
! SET SQ = (A(4)/S)**2 * (C**3/27 + D**2/4)
s = MAX(ABS(a(1)), ABS(a(2)), ABS(a(3)))
p1 = a(3)/(THREE*s)
q1 = a(2)/s
r1 = a(1)/s
t1 = q - 2.25D0*p*p
IF (ABS(t1) <= tol*ABS(q)) t1 = ZERO
w = FOURTH*r1*r1
w1 = HALF*p1*r1*t
w2 = q1*q1*t1/27.0D0
IF (w1 >= ZERO) THEN
w = w + w1
sq = w + w2
ELSE IF (w2 < ZERO) THEN
sq = w + (w1 + w2)
ELSE
w = w + w2
sq = w + w1
END IF
IF (ABS(sq) <= tol*w) sq = ZERO
rq = ABS(s/a(4))*SQRT(ABS(sq))
IF (sq >= ZERO) GO TO 40
! ALL ROOTS ARE REAL
arg = ATAN2(rq, -HALF*d)
cf = COS(arg/THREE)
sf = SIN(arg/THREE)
rt = SQRT(-c/THREE)
y1 = TWO*rt*cf
y2 = -rt*(cf + rt3*sf)
y3 = -(d/y1)/y2
x1 = y1 - p
x2 = y2 - p
x3 = y3 - p
IF (ABS(x1) > ABS(x2)) CALL Swap(x1,x2)
IF (ABS(x2) > ABS(x3)) CALL Swap(x2,x3)
IF (ABS(x1) > ABS(x2)) CALL Swap(x1,x2)
w = x3
IF (ABS(x2) < 0.1D0*ABS(x3)) GO TO 70
IF (ABS(x1) < 0.1D0*ABS(x2)) x1 = - (r/x3)/x2
z(1) = CMPLX(x1, ZERO,DP)
z(2) = CMPLX(x2, ZERO,DP)
z(3) = CMPLX(x3, ZERO,DP)
RETURN
! REAL AND COMPLEX ROOTS
40 ra =CubeRoot(-HALF*d - SIGN(rq,d))
rb = -c/(THREE*ra)
t = ra + rb
w = -p
x = -p
IF (ABS(t) <= tol*ABS(ra)) GO TO 41
w = t - p
x = -HALF*t - p
IF (ABS(x) <= tol*ABS(p)) x = ZERO
41 t = ABS(ra - rb)
y = HALF*rt3*t
IF (t <= tol*ABS(ra)) GO TO 60
IF (ABS(x) < ABS(y)) GO TO 50
s = ABS(x)
t = y/x
GO TO 51
50 s = ABS(y)
t = x/y
51 IF (s < 0.1D0*ABS(w)) GO TO 70
w1 = w/s
sum = ONE + t*t
IF (w1*w1 < 0.01D0*sum) w = - ((r/sum)/s)/s
z(1) = CMPLX(w, ZERO,DP)
z(2) = CMPLX(x, y,DP)
z(3) = CMPLX(x,-y,DP)
RETURN
! AT LEAST TWO ROOTS ARE EQUAL
60 IF (ABS(x) < ABS(w)) GO TO 61
IF (ABS(w) < 0.1D0*ABS(x)) w = - (r/x)/x
z(1) = CMPLX(w, ZERO,DP)
z(2) = CMPLX(x, ZERO,DP)
z(3) = z(2)
RETURN
61 IF (ABS(x) < 0.1D0*ABS(w)) GO TO 70
z(1) = CMPLX(x, ZERO,DP)
z(2) = z(1)
z(3) = CMPLX(w, ZERO,DP)
RETURN
! HERE W IS MUCH LARGER IN MAGNITUDE THAN THE OTHER ROOTS.
! AS A RESULT, THE OTHER ROOTS MAY BE EXCEEDINGLY INACCURATE
! BECAUSE OF ROUNDOFF ERROR. TO DEAL WITH THIS, A QUADRATIC
! IS FORMED WHOSE ROOTS ARE THE SAME AS THE SMALLER ROOTS OF
! THE CUBIC. THIS QUADRATIC IS THEN SOLVED.
! THIS CODE WAS WRITTEN BY WILLIAM L. DAVIS (NSWC).
70 aq(1) = a(1)
aq(2) = a(2) + a(1)/w
aq(3) = -a(4)*w
CALL QuadraticRoots(aq, z)
z(3) = CMPLX(w, ZERO,DP)
IF (AIMAG(z(1)) == ZERO) RETURN
z(3) = z(2)
z(2) = z(1)
z(1) = CMPLX(w, ZERO,DP)
RETURN
!-----------------------------------------------------------------------
! CASE WHEN D = 0
110 z(1) = CMPLX(-p, ZERO,DP)
w = SQRT(ABS(c))
IF (c < ZERO) GO TO 120
z(2) = CMPLX(-p, w,DP)
z(3) = CMPLX(-p,-w,DP)
RETURN
120 IF (p /= ZERO) GO TO 130
z(2) = CMPLX(w, ZERO,DP)
z(3) = CMPLX(-w, ZERO,DP)
RETURN
130 x = -(p + SIGN(w,p))
z(3) = CMPLX(x, ZERO,DP)
t = THREE*a(1)/(a(3)*x)
IF (ABS(p) > ABS(t)) GO TO 131
z(2) = CMPLX(t, ZERO,DP)
RETURN
131 z(2) = z(1)
z(1) = CMPLX(t, ZERO,DP)
RETURN
END Subroutine CubicRoots ! -----------------------------------------------
!+
SUBROUTINE QuarticRoots(a,z)
!----------------------------------------------------------------------------
! PURPOSE - Compute the roots of the real polynomial
! A(1) + A(2)*Z + ... + A(5)*Z**4
REAL, INTENT(IN) :: a(:)
COMPLEX, INTENT(OUT) :: z(:)
COMPLEX(DP) :: w
REAL(DP):: b,b2, c, d, e, h, p, q, r, t
REAL(DP),DIMENSION(4):: temp
REAL(DP):: u, v, v1, v2, x, x1, x2, x3, y
! NOTE - It is assumed that a(5) is non-zero. No test is made here
!----------------------------------------------------------------------------
IF (a(1)==0.0) THEN
z(1) = CZERO ! one root is obviously zero
CALL CubicRoots(a(2:), z(2:))
RETURN
END IF
b = a(4)/(FOUR*a(5))
c = a(3)/a(5)
d = a(2)/a(5)
e = a(1)/a(5)
b2 = b*b
p = HALF*(c - 6.0D0*b2)
q = d - TWO*b*(c - FOUR*b2)
r = b2*(c - THREE*b2) - b*d + e
! SOLVE THE RESOLVENT CUBIC EQUATION. THE CUBIC HAS AT LEAST ONE
! NONNEGATIVE REAL ROOT. IF W1, W2, W3 ARE THE ROOTS OF THE CUBIC
! THEN THE ROOTS OF THE ORIGINIAL EQUATION ARE
! Z = -B + CSQRT(W1) + CSQRT(W2) + CSQRT(W3)
! WHERE THE SIGNS OF THE SQUARE ROOTS ARE CHOSEN SO
! THAT CSQRT(W1) * CSQRT(W2) * CSQRT(W3) = -Q/8.
temp(1) = -q*q/64.0D0
temp(2) = 0.25D0*(p*p - r)
temp(3) = p
temp(4) = ONE
CALL CubicRoots(real(temp),z)
IF (AIMAG(z(2)) /= ZERO) GO TO 60
! THE RESOLVENT CUBIC HAS ONLY REAL ROOTS
! REORDER THE ROOTS IN INCREASING ORDER
x1 = DBLE(z(1))
x2 = DBLE(z(2))
x3 = DBLE(z(3))
IF (x1 > x2) CALL Swap(x1,x2)
IF (x2 > x3) CALL Swap(x2,x3)
IF (x1 > x2) CALL Swap(x1,x2)
u = ZERO
IF (x3 > ZERO) u = SQRT(x3)
IF (x2 <= ZERO) GO TO 41
IF (x1 >= ZERO) GO TO 30
IF (ABS(x1) > x2) GO TO 40
x1 = ZERO
30 x1 = SQRT(x1)
x2 = SQRT(x2)
IF (q > ZERO) x1 = -x1
temp(1) = (( x1 + x2) + u) - b
temp(2) = ((-x1 - x2) + u) - b
temp(3) = (( x1 - x2) - u) - b
temp(4) = ((-x1 + x2) - u) - b
CALL SelectSort(temp)
IF (ABS(temp(1)) >= 0.1D0*ABS(temp(4))) GO TO 31
t = temp(2)*temp(3)*temp(4)
IF (t /= ZERO) temp(1) = e/t
31 z(1) = CMPLX(temp(1), ZERO,DP)
z(2) = CMPLX(temp(2), ZERO,DP)
z(3) = CMPLX(temp(3), ZERO,DP)
z(4) = CMPLX(temp(4), ZERO,DP)
RETURN
40 v1 = SQRT(ABS(x1))
v2 = ZERO
GO TO 50
41 v1 = SQRT(ABS(x1))
v2 = SQRT(ABS(x2))
IF (q < ZERO) u = -u
50 x = -u - b
y = v1 - v2
z(1) = CMPLX(x, y,DP)
z(2) = CMPLX(x,-y,DP)
x = u - b
y = v1 + v2
z(3) = CMPLX(x, y,DP)
z(4) = CMPLX(x,-y,DP)
RETURN
! THE RESOLVENT CUBIC HAS COMPLEX ROOTS
60 t = DBLE(z(1))
x = ZERO
IF (t < ZERO) THEN
GO TO 61
ELSE IF (t == ZERO) THEN
GO TO 70
ELSE
GO TO 62
END IF
61 h = ABS(DBLE(z(2))) + ABS(AIMAG(z(2)))
IF (ABS(t) <= h) GO TO 70
GO TO 80
62 x = SQRT(t)
IF (q > ZERO) x = -x
70 w = SQRT(z(2))
u = TWO*DBLE(w)
v = TWO*ABS(AIMAG(w))
t = x - b
x1 = t + u
x2 = t - u
IF (ABS(x1) <= ABS(x2)) GO TO 71
t = x1
x1 = x2
x2 = t
71 u = -x - b
h = u*u + v*v
IF (x1*x1 < 0.01D0*MIN(x2*x2,h)) x1 = e/(x2*h)
z(1) = CMPLX(x1, ZERO,DP)
z(2) = CMPLX(x2, ZERO,DP)
z(3) = CMPLX(u, v,DP)
z(4) = CMPLX(u,-v,DP)
RETURN
80 v = SQRT(ABS(t))
z(1) = CMPLX(-b, v,DP)
z(2) = CMPLX(-b,-v,DP)
z(3) = z(1)
z(4) = z(2)
RETURN
END Subroutine QuarticRoots
!+
SUBROUTINE SelectSort(a)
! ---------------------------------------------------------------------------
! PURPOSE - Reorder the elements of in increasing order.
REAL(DP),INTENT(IN OUT),DIMENSION(:):: a
INTEGER:: j
INTEGER,DIMENSION(1):: k
! NOTE - This is a n**2 method. It should only be used for small arrays. <25
!----------------------------------------------------------------------------
DO j=1,SIZE(a)-1
k=MINLOC(a(j:))
IF (j /= k(1)) CALL Swap(a(k(1)),a(j))
END DO
RETURN
END Subroutine SelectSort ! -----------------------------------------------
!+
SUBROUTINE SolvePolynomial(quarticCoeff, cubicCoeff, quadraticCoeff, &
linearCoeff, constantCoeff, code, root1,root2,root3,root4)
! ---------------------------------------------------------------------------
REAL,INTENT(IN):: quarticCoeff
REAL,INTENT(IN):: cubicCoeff, quadraticCoeff
REAL,INTENT(IN):: linearCoeff, constantCoeff
INTEGER,INTENT(OUT):: code
COMPLEX,INTENT(OUT):: root1,root2,root3,root4
REAL,DIMENSION(5):: a
COMPLEX,DIMENSION(5):: z
!----------------------------------------------------------------------------
a(1)=constantCoeff
a(2)=linearCoeff
a(3)=quadraticCoeff
a(4)=cubicCoeff
a(5)=quarticCoeff
IF (quarticCoeff /= 0.) THEN
CALL QuarticRoots(a,z)
ELSE IF (cubicCoeff /= 0.) THEN
CALL CubicRoots(a,z)
ELSE IF (quadraticCoeff /= 0.) THEN
CALL QuadraticRoots(a,z)
ELSE IF (linearCoeff /= 0.) THEN
z(1)=CMPLX(-constantCoeff/linearCoeff, 0, DP)
outputCode=1
ELSE
outputCode=0 ! { no roots }
END IF
code=outputCode
IF (outputCode > 0) root1=z(1)
IF (outputCode > 1) root2=z(2)
IF (outputCode > 23) root3=z(3)
IF (outputCode > 99) root4=z(4)
RETURN
END Subroutine SolvePolynomial ! ------------------------------------------
SUBROUTINE SwapDouble(a,b)
! ---------------------------------------------------------------------------
! PURPOSE - Interchange the contents of a and b
REAL(DP),INTENT(IN OUT):: a,b
REAL(DP):: t
!----------------------------------------------------------------------------
t=b
b=a
a=t
RETURN
END Subroutine SwapDouble ! -----------------------------------------------
SUBROUTINE SwapSingle(a,b)
! ---------------------------------------------------------------------------
! PURPOSE - Interchange the contents of a and b
REAL(SP),INTENT(IN OUT):: a,b
REAL(SP):: t
!----------------------------------------------------------------------------
t=b
b=a
a=t
RETURN
END Subroutine SwapSingle ! -----------------------------------------------
END Module PolynomialRoots ! ==============================================