!###########################################################
! page.79-80 (例題4)
! 1次元ニュートン・ラフソン法の複素変数プログラム
!
!###########################################################
PROGRAM prog16
IMPLICIT REAL(8)(a-h,o-z)
IMPLICIT INTEGER(i-n)

INTEGER,PARAMETER :: KMAX=10
COMPLEX :: z,fz,zk

!---------------
z=CMPLX(-1.d0,-3.d0)
delta=1.d-8
!---------------

fz=funcz(z)
k=0
WRITE(*,"('k=',i2,SP,', z=',1P2E14.6,'*i, f(z)=',1P2E14.6,'*i')") k,REAL(z),AIMAG(z),REAL(fz),AIMAG(fz)

DO k=1,KMAX

zk=fnewton(z)
fz=funcz(zk)
WRITE(*,"('k=',i2,SP,', z=',1P2E14.6,'*i, f(z)=',1P2E14.6,'*i')") k,REAL(zk),AIMAG(zk),REAL(fz),AIMAG(fz)

IF(cabs(fz)**2<delta)GOTO 100

z=zk

END DO
WRITE(*,*) '解が求まりませんでした'
STOP

100 CONTINUE
WRITE(*,*)
WRITE(*,"('解は z=',SP,1P2E14.6,'*i')") REAL(zk),AIMAG(zk)
STOP

!###########################################################
CONTAINS

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
FUNCTION funcz(z)
   COMPLEX,INTENT(IN) :: z
   COMPLEX :: funcz
   funcz=z**3+6.d0*z**2+21.d0*z+32.d0
RETURN
END FUNCTION funcz
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  
FUNCTION fnewton(xx)
   COMPLEX,INTENT(IN) :: xx
   COMPLEX :: fnewton
fnewton=xx-fz/(3.d0*xx**2+12.d0*xx+21.d0)
RETURN
END FUNCTION fnewton
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!  

END PROGRAM prog16