SUBROUTINE POL_direct(N, coef, coef_1, X0, EPS, ITMAX, RAD, B, IER, n_div)
! Calculeeaza 1 radacina a lui p(x)=a_0+a_1*x+...+a_n*x**n
! Prin iterare in polinomul original
! Parametri:
! n: gradul polinomului
! A: vectorul coeficientilor: a0,a1,...,an
! x0: aproximatia initiala
! eps: toleranta pentru metoda Newton
! itmax: -INTRARE: Numarul maxim de iteratii
!		 -IESIRE:  Numarul efectiv de iteratii
! rad: radacina/ultima aproximatie
! B = p(x); C =q(x)
! ier: cod de intrare/iesire:
!		- INTRARE: tiparirea iteratiilor
!		- IESIRE: cod de incheiare a iteratiei:
!			0: toleranta 'eps' atinsa
!			1: numar iteratii depasit
!		    2: proces stationar (cu "period")
!			-1: divergenta
!			-2:	numitor nul.
! n_div: numar de divergente locale (pana la ultima iteratie)

!real A(0:N)
implicit real(8) (a-h,o-z)

real(8):: coef(0:n), coef_1(0:n-1)

logical stationar
integer width
logical print_it

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

logical print_div


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

print_it = ier /=0

print_div =.false.

! Antet iteratii:
if(print_it) then
write(1,*) 'Iteratia		x			Diferenta curenta'
endif

iter=1
dif_ant =100.
aa(0) =x0
n_div =0

200	Z=X0
b =polval(n, coef,z)
c =polval(n-1, coef_1,z)

! Test numitor nul:
!if(c ==0.0) then
if(abs(c) <1E-10) then

write(1,*); write(1,*) ' Radacina', z, 'Numitor:', c;
ier =-2

rad =z
itmax =iter
RETURN
endif


X1 =X0 -b/c

dif =x1 -x0

if(print_it) then
!write(1,*) iter, x1, dif
write(1, '(/, I9, t13, 1pG14.7, 1pG11.4, \)') iter, x1, dif
endif

! Test toleranta:
	IF(ABS(dif) <= EPS) THEN
	IER=0
	GO TO 500
	ENDIF

! Test proces stationar:
aa(iter) =x1		! 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', x1
	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', x1
	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

ier =2
goto 500
ENDIF

! Divergenta locala:
IF(abs(dif) > abs(dif_ant)) THEN
!ier =-2 
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 ==ITMAX) THEN
	IER=1
	GO TO 500

	ELSE
	iter =iter+1
	X0=X1
	dif_ant =dif
	GO TO 200
	ENDIF

500	RAD=X1
	itmax=iter

END
