subroutine Halley(SUB, x0, eps, lnit, rad, kod, n_div)
use work
! Solutia ecuatiei F(x)=0, prin metoda Halley.
! Parametri:
!	SUB	:numele subrutinei care calculeaza valoarea lui F si a derivatei F1.
!		 Se declara EXTERNAL in programul care apeleaza NEWTON.
!	X0	:aproximatia initiala
!	eps	:toleranta
!	lnit: -Intrare =nr.maxim de iteratii; -Iesire =nr.efectiv de iteratii.
!	rad = solutia ecuatiei
!	kod	: -INTRARE: cod de tiparire a iteratiilor: 0 = Nu; #0 = Da.
!       : -IESIRE : cod de terminare/eroare:
!		 0:	toleranta "eps" atinsa;
!		 1:	depasirea nr.maxim de iteratii admis "lnit";
!        2: Proces stationar;
!		-1: divergenta
!       -2:	derivata functiei in X0 este ~ nula.
!
!* Exemplu de subrutina SUB:
! subroutine SUB(f,f1,x)
! common p	! optional, in cazul cand F include un parametru "p".
!			! "p" se citeste in programul apelant.
! f =...
! f1=...
! end

logical print_it
logical stationar
integer width

real aa(0:lnit)	! Tablou pentru toate val. - pt. proces stationar
integer period, per

logical print_div
real num    ! numitorul in Formula Halley

!epsX(x) =1./2**(24 -exponent(x))

print_it = kod /=0

print_div =.false.

if(print_it) then
print *
write(1, '(a)') ' Iteratia    x'
endif


! Initialize:
aa =0.
aa(0) =x0

iter =0

x =x0
dif_ant =100.
!n_div =0   ! in Main

!eps_num =1E-12

100	iter =iter +1

if(not_Pol) then
call sub(f, f1, f2, x)

else
f  =polval(n, coef, x)
f1 =polval(n-1, coef_1, x)
f2 =polval(n-2, coef_2, x)
endif

num =f1*f1 -f2*f/2.

! Test numitor nul:
if(abs(num) < eps_num) then		! tiny(f1)

print *, 'Numitor =', num

write(1,*) 'Numitor =', num

kod=-2
goto 500
endif

!x  =x0 -f/f1		!dif
x =x0 -f1*f/(f1*f1 -f2*f/2.)    ! num

dif =x-x0

iwi =width(iter)
	if(print_it) then
	print 90, iter,x
	write(1,90) iter, x
	endif

90	format(/, I9, t13, 1pG14.7, \)

! EPS = atins:
!IF(abs(dif) <= eps) THEN
IF(abs(dif) <= eps) THEN
kod =0
goto 500
ENDIF

! Test Stationar:
aa(iter) =x		! pt. stationar
per =period(iter, aa, i)

IF(per >0) then
iwiter =width(iter)
iwi =width(i)
iwper =width(per)

!	write(*, '(/, a, 1pG14.7, \)') ' Radacina', x
	print '(/, a, I<iwiter>, a, I<iwi>, a)', ' x(', iter, ') = x(', i, ')'
	print '(a, I<iwper>)', ' ** Proces stationar de perioada ', per

!	write(1, '(//, a, 1pG14.7, \)') ' Radacina', x
	write(1, '(/, a, I<iwiter>, a1, \)') ' Iteratia ', iter,': '
	write(1, '(a, I<iwiter>, a, I<iwi>, a)') ' x(', iter, ') = x(', i, ')'
	write(1, '(a, I<iwper>)') ' ** Proces stationar de perioada = ', per

kod =2
goto 500
ENDIF

! Divergenta locala:
IF(abs(dif) > abs(dif_ant)) THEN
!kod =-1 
n_div =n_div+1

if(print_it) then
if(.not. print_div) print *, ' ** Divergenta locala la iteratia:'
print_div =.true.

print *, iter

write(1, '(a)') ' ** Divergenta locala'
endif
ENDIF

! Test numar de iteratii:
IF(iter < lnit) THEN
goto 200

ELSE
	if(abs(dif) > abs(dif_ant)) then
	kod =-1		! Divergenta
	else
	kod=1		! LNIT depasit
	endif

goto 500
ENDIF

! Noua iteratie
200 continue
dif_ant =dif
x0 =x
goto 100


500 continue
rad =x
lnit =iter

if(.not. print_it .and. kod ==2 .and. per <=50) then
write(1,*)
do ii =i, iter
write(1, '(I6, t14, G16.7)') ii, aa(ii)
enddo
endif

print '(/, a, 1pg14.7)','  Diferenta curenta = ', dif

if(print_it .or. kod ==2) write(1,'(/)')
write(1, '(a, 1pg14.7)') '  Diferenta curenta =', dif

END
