SUBROUTINE POL_direct(N,A,X0,EPS,ITMAX,RAD,B,IER)
! Calculeeaza 1 radacina a lui p(x)=a0+a1*x+...+an*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
!		   -1: proces stationar (x1=x0)
!		   -2: numitor nul

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

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

! antet iteratii:
if(ier /=0) then
write(1,*) '	Iteratia	x					Diferenta curenta'
endif

200	Z=X0
	C=A(N)
	B =A(N)

DO K=N-1,0,-1
b =a(k) +z*b
if(k ==0) exit
c =b +z*c
ENDDO

! Test numitor nul:
if(c ==0.0) then
ier =-2
rad =z
return
endif


X1 =X0 -b/c

dif =x1 -x0

if(ier /=0) then
write(1,*) iter, x1, dif
endif

! Test toleranta si proces stationar:
	IF(ABS(dif) <=EPS) THEN
	IER=0
	GO TO 500
	ELSE IF	(abs(dif) < ulp(x0)) then
	! Test proces stationar: x1 = x0
	ier =-1
	print *, 'x0=', x0,':   MaxEps =', ulp(x0),'  Diferenta =', dif
	write(1,*) ' x0=', x0,':   MaxEps =', ulp(x0),'  Diferenta =', dif
	goto 500
	ENDIF

! Test numar de itaratii
	IF(iter.EQ.ITMAX) THEN
	IER=1
	GO TO 500
	ELSE
	iter=iter+1
	X0=X1
	GO TO 200
	ENDIF

500	RAD=X1
	itmax=iter

END
