Fortran解一元四次方程程序

 一元四次方程， 沈天珩求根公式,  具体公式 详见百科链接；

程序由main.f90和输入文件abcde.para构成，abcde分别为

 程序的结果将给出abcde取值，参数符合天珩求根公式第几定理，四个根，以及将四个根分别带回方程验证。在下面会给出示例。

1.main.f90

program main
    IMPLICIT NONE
integer :: i1
real*8, parameter :: pi = acos(-1.0d0)
complex*16 :: X(4), test(4)
complex*16 :: X1,X2,X3,X4
real*8 :: a,b,c,d,e
real*8 :: DD ,EE, FF, AA, BB, CC, delta
real*8 :: z1, z2, z
real*8 :: y1, y2, y3
real*8 :: g3
complex*16,parameter :: i=(0.0d0,1.0d0)
real*8 :: theta
open(10,file="abcde.para")
read(10,*)a
read(10,*)b
read(10,*)c
read(10,*)d
read(10,*)e
close(10)
write(*,*)"The a is",a,"The b is",b
write(*,*)"The c is",c,"The d is",d
write(*,*)"The e is",e

DD=3.0d0*b**2.0d0 - 8.0d0*a*c

EE=-b**3.0d0 + 4.0d0*a*b*c - 8.0d0*a**2.0d0*d

FF=3.0d0*b**4.0d0 + 16.0d0*a**2.0d0*c**2.0d0 - 16.0d0*a*b**2.0d0*c + 16.0d0*a**2.0d0*b*d - 64.0d0*a**3.0d0*e

AA=DD**2.0d0 - 3.0d0*FF

BB=DD*FF - 9.0d0*EE**2.0d0

CC=FF**2.0d0 - 3.0d0*DD*EE**2.0d0

delta=BB**2.0d0 - 4.0d0*AA*CC

if (FF==0.0 .AND. EE==0.0 .AND. DD==0.0) THEN
WRITE(*,*)"Tianheng formula categories NO1"

X1=-4.0*e/d
X2=-4.0*e/d
X3=-4.0*e/d
X4=-4.0*e/d

else if (DD*EE*FF /=0.0 .and. CC==0.0 .AND. BB==0.0 .AND. AA==0.0) THEN
WRITE(*,*)"Tianheng formula categories NO2"
X1 = (-b*DD+9*EE)/(4*a*DD)
X2 = (-b*DD-3*EE)/(4*a*DD)
X3 = (-b*DD-3*EE)/(4*a*DD)
X4 = (-b*DD-3*EE)/(4*a*DD)

else if (FF==0.0 .AND. EE==0.0 .AND. DD>0.0)THEN
WRITE(*,*)"Tianheng formula categories NO3"
X1 = (-b+sqrt(DD))/(4*a)
X2 = (-b+sqrt(DD))/(4*a)
X3 = (-b-sqrt(DD))/(4*a)
X4 = (-b-sqrt(DD))/(4*a)

else if (FF==0.0 .AND. EE==0.0 .AND. DD<0.0)THEN
WRITE(*,*)"Tianheng formula categories NO4"
X1 = (-b+i*sqrt(-DD))/(4*a)
X2 = (-b+i*sqrt(-DD))/(4*a)
X3 = (-b-i*sqrt(-DD))/(4*a)
X4 = (-b-i*sqrt(-DD))/(4*a)

else if (AA*BB*CC /=0.0 .AND. delta==0.0 .and. AA*BB>0.0) THEN
WRITE(*,*)"Tianheng formula categories NO5 bingo"
X1=(-b+2*AA*EE/BB+SQRT(2.0*BB/AA))/(4*a)
X2=(-b+2*AA*EE/BB-SQRT(2.0*BB/AA))/(4*a)
X3=(-b-2*AA*EE/BB)/(4*a)
X4=(-b-2*AA*EE/BB)/(4*a)

else if (AA*BB*CC /=0.0 .AND. delta==0.0 .and. AA*BB<0.0) THEN
WRITE(*,*)"Tianheng formula categories NO6"
X1=(-b+2*AA*EE/BB+i*SQRT(-2.0*BB/AA))/(4*a)
X2=(-b+2*AA*EE/BB-i*SQRT(-2.0*BB/AA))/(4*a)
X3=(-b-2*AA*EE/BB)/(4*a)
X4=(-b-2*AA*EE/BB)/(4*a)


else if (delta>0.0) THEN
WRITE(*,*)"Tianheng formula categories NO7 bingo"
z1=AA*DD + 3.0d0*(-BB+sqrt(delta))/2.0d0
z2=AA*DD + 3.0d0*(-BB-sqrt(delta))/2.0d0
z=DD**2.0d0 - DD*(g3(z1)+g3(z2)) + (g3(z1)+g3(z2))**2.0d0 -3.0d0*AA
 
X1=(-b + sign(1.0d0,EE) * ((DD + g3(z1) + g3(z2))/3.0d0)**(1.0d0/2.0d0) &
 + ((2.0d0*DD - g3(z1) - g3(z2) + 2.0d0*z**(1.0d0/2.0d0))/3.0d0)**(1.0d0/2.0d0))/(4.0d0*a)

X2=(-b + sign(1.0d0,EE) * ((DD + g3(z1) + g3(z2))/3.0d0)**(1.0d0/2.0d0) &
 - ((2.0d0*DD - g3(z1) - g3(z2) + 2.0d0*z**(1.0d0/2.0d0))/3.0d0)**(1.0d0/2.0d0))/(4.0d0*a)

X3=(-b - sign(1.0d0,EE) * ((DD + g3(z1) + g3(z2))/3.0d0)**(1.0d0/2.0d0) & 
 + i* real((-2.0d0*DD + g3(z1) + g3(z2) + 2.0d0*z**(1.0d0/2.0d0))/3.0d0)**(1.0d0/2.0d0))/(4.0d0*a)

X4=(-b - sign(1.0d0,EE) * ((DD + g3(z1) + g3(z2))/3.0d0)**(1.0d0/2.0d0) &
 - i* real((-2.0d0*DD + g3(z1) + g3(z2) + 2.0d0*z**(1.0d0/2.0d0))/3.0d0)**(1.0d0/2.0d0))/(4.0d0*a)

else if (delta<0.0) THEN
	theta = dacos( (3.0d0*BB - 2.0d0*AA*DD)/(2.0d0*AA*sqrt(AA)) ) 

	y1=( DD-2.0d0*sqrt(AA)*dcos(theta/3.0d0) )/3.0d0

	y2=( DD+sqrt(AA)*(dcos(theta/3.0d0)+ SQRT(3.0D0)*dsin(theta/3.0d0) ))/3.0d0

	y3=( DD+sqrt(AA)*(dcos(theta/3.0d0)- SQRT(3.0D0)*dsin(theta/3.0d0) ))/3.0d0



 	if (EE==0.0d0 .and. DD>0.0d0 .and. FF>0.0d0) then
	WRITE(*,*)"Tianheng formula categories NO8<1>"
	X1=(-b + sqrt( DD+ 2.0d0*sqrt(FF)) )/(4.0d0*a)
	X2=(-b - sqrt( DD+ 2.0d0*sqrt(FF)) )/(4.0d0*a)
	X3=(-b + sqrt( DD- 2.0d0*sqrt(FF)) )/(4.0d0*a)
	X4=(-b - sqrt( DD- 2.0d0*sqrt(FF)) )/(4.0d0*a)

 	elseif (EE==0.0d0 .and. DD<0.0d0 .and. FF>0.0d0) then
	WRITE(*,*)"Tianheng formula categories NO8<2>"
	X1=(-b + sqrt( -DD+ 2.0d0*sqrt(FF)) )/(4.0d0*a)
	X2=(-b - sqrt( -DD+ 2.0d0*sqrt(FF)) )/(4.0d0*a)
	X3=-b/(4.0d0*a) + i*sqrt( -DD- 2.0d0*sqrt(FF)) /(4.0d0*a)
	X4=-b/(4.0d0*a) - i*sqrt( -DD- 2.0d0*sqrt(FF)) /(4.0d0*a)

 	elseif (EE==0.0d0 .and. FF<0.0d0) then
	WRITE(*,*)"Tianheng formula categories NO8<3>"
	X1=(-2.0d0*b + sqrt(2.0d0*DD + 2.0d0*sqrt(AA-FF)))/(8.0d0*a) + i*sqrt(-2.0d0*DD+2.0d0*sqrt(AA-FF))/(8.0d0*a)
	X2=(-2.0d0*b + sqrt(2.0d0*DD + 2.0d0*sqrt(AA-FF)))/(8.0d0*a) - i*sqrt(-2.0d0*DD+2.0d0*sqrt(AA-FF))/(8.0d0*a)
	X3=(-2.0d0*b - sqrt(2.0d0*DD + 2.0d0*sqrt(AA-FF)))/(8.0d0*a) + i*sqrt(-2.0d0*DD+2.0d0*sqrt(AA-FF))/(8.0d0*a)
	X4=(-2.0d0*b - sqrt(2.0d0*DD + 2.0d0*sqrt(AA-FF)))/(8.0d0*a) - i*sqrt(-2.0d0*DD+2.0d0*sqrt(AA-FF))/(8.0d0*a)

 	elseif (EE/=0.0d0 .and. DD>0.0d0 .and. FF>0.0d0) then
	WRITE(*,*)"Tianheng formula categories NO8<4>(1)"
	X1=(-b +  sign(1.0d0,EE)*sqrt(y1) + (sqrt(y2) + sqrt(y3)) )/(4.0d0*a)
	X2=(-b +  sign(1.0d0,EE)*sqrt(y1) - (sqrt(y2)+sqrt(y3)) )/(4.0d0*a)
	X3=(-b -  sign(1.0d0,EE)*sqrt(y1) + (sqrt(y2)-sqrt(y3)) )/(4.0d0*a)
	X4=(-b -  sign(1.0d0,EE)*sqrt(y1) - (sqrt(y2)-sqrt(y3)) )/(4.0d0*a)

 	elseif (EE/=0.0d0 .and. DD<0.0d0 .or. FF<0.0d0) then
	WRITE(*,*)"Tianheng formula categories NO8<4>(2)"
	X1=(-b -sqrt(y2))/(4.0d0*a) + i*(sign(1.0d0,EE)*sqrt(-y1)+sqrt(-y3))/(4.0*a)
	X2=(-b -sqrt(y2))/(4.0d0*a) - i*(sign(1.0d0,EE)*sqrt(-y1)+sqrt(-y3))/(4.0*a)
	X3=(-b +sqrt(y2))/(4.0d0*a) + i*(sign(1.0d0,EE)*sqrt(-y1)-sqrt(-y3))/(4.0*a)
	X4=(-b +sqrt(y2))/(4.0d0*a) - i*(sign(1.0d0,EE)*sqrt(-y1)-sqrt(-y3))/(4.0*a)
	endif
 

ENDIF

X(1)=x1
X(2)=X2
X(3)=X3
X(4)=X4

do i1=1,4
   test(i1)= a*X(i1)**4.0d0 + b*X(i1)**3.0d0 + c*X(i1)**2.0d0 + d*X(i1)**1.0d0 + e
enddo


WRITE(*,*)" "
WRITE(*,*)"the X1 is", real(X(1)) ,"+", imag(X(1)) ,"i"
WRITE(*,*)"the X2 is", real(X(2)) ,"+", imag(X(2)) ,"i"
WRITE(*,*)"the X3 is", real(X(3)) ,"+", imag(X(3)) ,"i"
WRITE(*,*)"the X4 is", real(X(4)) ,"+", imag(X(4)) ,"i"
WRITE(*,*)" "
WRITE(*,*)"Bringing x1 back to the original equation tests as",real(test(1))
WRITE(*,*)"Bringing x2 back to the original equation tests as",real(test(2))
WRITE(*,*)"Bringing x3 back to the original equation tests as",real(test(3))
WRITE(*,*)"Bringing x4 back to the original equation tests as",real(test(4))
WRITE(*,*)" "

END

function g3(number1)
real*8 :: number1
real*8 :: g3

g3=sign(1.0,real(number1))*abs(number1)**(1.0d0/3.0d0)
return
end

2.abcde.para

1.0	!a
2.0	!b
2.0	!c
30.0	!d
50.0	!e

3.result