C **********************************************************25.03.98 C Test program for QUINTIC (solving quartic equations) C INPUT: iTestMethod 0 : input of coefficients C 1 : input of four real roots C NO WARRANTY, ALWAYS TEST THIS PROGRAM AFTER DOWNLOADING C ****************************************************************** program QuinticTest implicit double precision (a-h,o-z) dimension dd(0:5),sol(5),soli(5) common /start/ xstart C C input roots directly f = (a-x)*(b-x)*(c-x)*(d-x)*(e-x) read(5,*)iTestMethod, xstart if (iTestMethod.eq.1) then read(5,*) a,b,c,d,e dd(5) = -1.d+0 dd(4) = a + b + c + d + e dd(3) = -( a*b + (a+b)*c + (a+b+c)*d + (a+b+c+d)*e ) dd(2) = a*b*c + (a*b+(a+b)*c)*d + (a*b+(a+b)*c+(a+b+c)*d)*e dd(1) = -( a*b*c*d + (a*b*c + (a*b + (a+b)*c)*d)*e ) dd(0) = a*b*c*d*e write(6,*)'QUINTIC test:' write(6,*)'input of four real roots:' write(6,'(a5,f16.8)')'x1 = ',a write(6,'(a5,f16.8)')'x2 = ',b write(6,'(a5,f16.8)')'x3 = ',c write(6,'(a5,f16.8)')'x4 = ',d write(6,'(a5,f16.8)')'x5 = ',e else do 10 i=0,5 10 read(5,*) dd(i) write(6,*)'QUINTIC test:' write(6,*)'input of coefficients:' endif do 15 i=0,5 15 write(6,'(i5,f16.8)')i,dd(i) C call quintic(dd,sol,soli,Nsol) write(6,*)'solutions:' write(6,'(a5,2a16)')' #','real','imaginary' do 20 i=1,5 20 write(*,'(i5,2f16.8)') i,sol(i),soli(i) C C check real solutions write(6,*)'check real solution:' write(6,'(a5,2a16)')' #','root','f(root)' do 30 i=1,Nsol y = dd(0) do 31 k=1,5 31 y = y + dd(k)*sol(i)**k write(6,'(i5,2f16.8)')i,sol(i),y 30 continue C open(8,file='quintic.chk',status='unknown',access='sequential') do 40 i=0,5 40 write(8,'(a2,f16.8,a7)') ' (',dd(i),', 0.0 )' do 50 i=1,5 50 write(8,'(a2,f16.8,a1,f16.8,a1)') ' (',sol(i),',',soli(i),')' close(8) C stop end C ***QUINTIC************************************************25.03.98 C Solution of a quintic equation by a hybrid method: C first real solution is obtained numerically by the Newton method, C the remaining four roots are obtained analytically by QUARTIC C NO WARRANTY, ALWAYS TEST THIS SUBROUTINE AFTER DOWNLOADING C ****************************************************************** C dd(0:4) (i) vector containing the polynomial coefficients C sol(1:4) (o) results, real part C soli(1:4) (o) results, imaginary part C Nsol (o) number of real solutions C ================================================================== subroutine quintic(dd,sol,soli,Nsol) implicit double precision (a-h,o-z) parameter( eps = 1.d-8 ) dimension dd(0:5),sol(5),soli(5) dimension dd4(0:4),sol4(4),soli4(4) common /start/ xstart C Nsol = 0 if (dd(5).eq.0.d+0) then write(6,*)'ERROR: NOT A QUINTIC EQUATION' return endif C C Newton iteration of one real root xs = xstart 4 sum = dd(0) do 5 i=1,5 5 sum = sum + dd(i)*xs**i sum1 = dd(1) do 6 i=1,4 6 sum1 = sum1 + float(i+1)*dd(i+1)*xs**i xnew = xs - sum/sum1 if (abs(xnew-xs).gt.eps) then xs = xnew goto 4 endif e = xnew C C "e" is one real root of quintic equation C reduce quintic to quartic equation using "e" dd4(4) = dd(5) do 8 i=4,1,-1 8 dd4(i-1) = dd(i) + e*dd4(i) do 9 i=0,4 9 write(6,*) dd4(i) call quartic(dd4,sol4,soli4,Nsol4) sol(1) = e soli(1) = 0.d+0 do 10 i=1,4 sol(i+1) = sol4(i) soli(i+1) = soli4(i) 10 continue Nsol = Nsol4 + 1 C return END C ***QUARTIC************************************************25.03.98 C Solution of a quartic equation C ref.: J. E. Hacke, Amer. Math. Monthly, Vol. 48, 327-328, (1941) C NO WARRANTY, ALWAYS TEST THIS SUBROUTINE AFTER DOWNLOADING C ****************************************************************** C dd(0:4) (i) vector containing the polynomial coefficients C sol(1:4) (o) results, real part C soli(1:4) (o) results, imaginary part C Nsol (o) number of real solutions C ================================================================== subroutine quartic(dd,sol,soli,Nsol) implicit double precision (a-h,o-z) dimension dd(0:4),sol(4),soli(4) dimension AA(0:3),z(3) C Nsol = 0 if (dd(4).eq.0.d+0) then write(6,*)'ERROR: NOT A QUARTIC EQUATION' return endif a = dd(4) b = dd(3) c = dd(2) d = dd(1) e = dd(0) C p = (-3.d+0*b**2 + 8.d+0*a*c)/(8.d+0*a**2) q = (b**3 - 4.d+0*a*b*c + 8.d+0*d*a**2)/(8.d+0*a**3) r = (-3.d+0*b**4 + 16.d+0*a*b**2*c - 64.d+0*a**2*b*d + & 256.d+0*a**3*e)/(256.d+0*a**4) C C solve cubic resolvent AA(3) = 8.d+0 AA(2) = -4.d+0*p AA(1) = -8.d+0*r AA(0) = 4.d+0*p*r - q**2 call cubic(AA,z,ncube) C zsol = -1.d+99 do 5 i=1,ncube 5 zsol = max(zsol,z(i)) z(1) = zsol xK2 = 2.d+0 * z(1) - p xK = sqrt(xK2) xL = q/(2.d+0 * xK) sqp = xK2 - 4.d+0*(z(1) + xL) sqm = xK2 - 4.d+0*(z(1) - xL) C do 10 i=1,4 10 soli(i) = 0.d+0 if (sqp.ge.0.d+0 .and. sqm.ge.0.d+0) then sol(1) = 0.5d+0*( xK + sqrt(sqp)) sol(2) = 0.5d+0*( xK - sqrt(sqp)) sol(3) = 0.5d+0*(-xK + sqrt(sqm)) sol(4) = 0.5d+0*(-xK - sqrt(sqm)) Nsol = 4 else if (sqp.ge.0.d+0 .and. sqm.lt.0.d+0) then sol(1) = 0.5d+0*(xK + sqrt(sqp)) sol(2) = 0.5d+0*(xK - sqrt(sqp)) sol(3) = -0.5d+0*xK sol(4) = -0.5d+0*xK soli(3) = sqrt(-0.25d+0 * sqm) soli(4) = -sqrt(-0.25d+0 * sqm) Nsol = 2 else if (sqp.lt.0.d+0 .and. sqm.ge.0.d+0) then sol(1) = 0.5d+0*(-xK + sqrt(sqm)) sol(2) = 0.5d+0*(-xK - sqrt(sqm)) sol(3) = 0.5d+0*xK sol(4) = 0.5d+0*xK soli(3) = sqrt(-0.25d+0 * sqp) soli(4) = -sqrt(-0.25d+0 * sqp) Nsol = 2 else if (sqp.lt.0.d+0 .and. sqm.lt.0.d+0) then sol(1) = -0.5d+0*xK sol(2) = -0.5d+0*xK soli(1) = sqrt(-0.25d+0 * sqm) soli(2) = -sqrt(-0.25d+0 * sqm) sol(3) = 0.5d+0*xK sol(4) = 0.5d+0*xK soli(3) = sqrt(-0.25d+0 * sqp) soli(4) = -sqrt(-0.25d+0 * sqp) Nsol = 0 endif do 20 i=1,4 20 sol(i) = sol(i) - b/(4.d+0*a) C return END C ***CUBIC************************************************08.11.1986 C Solution of a cubic equation C Equations of lesser degree are solved by the appropriate formulas. C The solutions are arranged in ascending order. C NO WARRANTY, ALWAYS TEST THIS SUBROUTINE AFTER DOWNLOADING C ****************************************************************** C A(0:3) (i) vector containing the polynomial coefficients C X(1:L) (o) results C L (o) number of valid solutions (beginning with X(1)) C ================================================================== SUBROUTINE CUBIC(A,X,L) IMPLICIT DOUBLE PRECISION(A-H,O-Z) DIMENSION A(0:3),X(3),U(3) PARAMETER(PI=3.1415926535897932D+0,THIRD=1.D+0/3.D+0) INTRINSIC MIN,MAX,ACOS C C define cubic root as statement function CBRT(Z)=SIGN(ABS(Z)**THIRD,Z) C C ==== determine the degree of the polynomial ==== C IF (A(3).NE.0.D+0) THEN C C cubic problem W=A(2)/A(3)*THIRD P=(A(1)/A(3)*THIRD-W**2)**3 Q=-.5D+0*(2.D+0*W**3-(A(1)*W-A(0))/A(3)) DIS=Q**2+P IF (DIS.LT.0.D+0) THEN C three real solutions! C Confine the argument of ACOS to the interval [-1;1]! PHI=ACOS(MIN(1.D+0,MAX(-1.D+0,Q/SQRT(-P)))) P=2.D+0*(-P)**(5.D-1*THIRD) DO 100 I=1,3 100 U(I)=P*COS((PHI+DBLE(2*I)*PI)*THIRD)-W X(1)=MIN(U(1),U(2),U(3)) X(2)=MAX(MIN(U(1),U(2)),MIN(U(1),U(3)),MIN(U(2),U(3))) X(3)=MAX(U(1),U(2),U(3)) L=3 ELSE C only one real solution! DIS=SQRT(DIS) X(1)=CBRT(Q+DIS)+CBRT(Q-DIS)-W L=1 END IF C ELSE IF (A(2).NE.0.D+0) THEN C C quadratic problem P=5.D-1*A(1)/A(2) DIS=P**2-A(0)/A(2) IF (DIS.GE.0.D+0) THEN C two real solutions! X(1)=-P-SQRT(DIS) X(2)=-P+SQRT(DIS) L=2 ELSE C no real solution! L=0 END IF C ELSE IF (A(1).NE.0.D+0) THEN C C linear equation X(1)=-A(0)/A(1) L=1 C ELSE C no equation L=0 END IF C C ==== perform one step of a newton iteration in order to minimize C round-off errors ==== DO 110 I=1,L X(I)=X(I)-(A(0)+X(I)*(A(1)+X(I)*(A(2)+X(I)*A(3)))) * /(A(1)+X(I)*(2.D+0*A(2)+X(I)*3.D+0*A(3))) 110 CONTINUE RETURN END