! u =f/f'
! om1 =f'
! om2 =f'(x+alfa*u)
! om3 =f'(x +beta*u +gama*f/om2)
!
! Ordin 3:
! j =f/(a1*om1+a2*om2)
! Ordin 4, 5:
! j =f/(a1*om1+a2*om2+a3*om3)
! (j --> jf)
!
! Metoda:
! x_(n+1) =x_n -j(x_n)
!
subroutine JARRATT(x0, eps, lnit, rad, kod, n_div)
use work
! Solutia ecuatiei F(x)=0, prin metoda JARRATT.
! Parametri:
!	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.
!
!  SUB: subrutina care calculeaza valoarea lui f, a derivatei f1, si a functiilor: om1=f1, om2, om3.
!* Exemplu de subrutina SUB:
! subroutine SUB(x, f, om1, om2, om3)
! use work
!
! f0(x0 =... (expresia lui f)
! f1(x) =... (expresia lui f')
! ...
! 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 Jarratt
real jf     ! functia de iterare "j"

!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(x, f, om1, om2, om3)

! om1 =f1
if(abs(om1) <eps_num) then
kod =-2
write(1,*) 'f'' =', om1
goto 500
endif

else
f  =polval(n, coef, x)
om1 =polval(n-1, coef_1, x)

if(abs(om1) <eps_num) then
kod =-2
write(1,*) 'f'' =', om1
goto 500
endif

u =f/om1
x2 =x +alfa*u
om2 =polval(n-1, coef_1, x2)

if(order >3) then
x3 =x +beta*u +gama*f/om2
om3 =polval(n-1, coef_1, x3)
endif

endif


!u =f/f'
! om1 =f'
! om2 =f'(x+alfa*u)
! om3 =f'(x +beta*u +gama*f/om2)

! Defineste "num":
! j =f/num
!

call gen_NUM(om1, om2, om3, num)

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

print *, 'Numitor =', num

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

kod=-2
goto 500
endif

jf =f/num
x =x0 -jf        ! u/num; u*factor;

dif =-jf     !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)

	print '(/, a, I<iwiter>, a, I<iwi>, a)', ' x(', iter, ') = x(', i, ')'
	print '(a, I<iwper>)', ' ** Proces stationar de perioada ', per

	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
