! y =x -f(x)/f'(x)
! z =y -f(x)/(f(x) -2*f(y))*f(y)/f'(x)
! w =z -f(z)*f[x,y]/div_dif(f,x,z)/div_dif(f,y,z)*(1+ f(z)/f(x))
! x1 =w -f(w)/( div_dif(f,x,w) +( div_dif_2(f,y,x,z) -div_dif_2(f,y,x,w) -div_dif_2(f,z,x,w) )*(x-w) )
!
! Ordin 15:

subroutine Pas_Ordin_15(x0, x, iter, print_pas, kod)
use work
! Solutia ecuatiei F(x)=0, prin metoda de Ordin 15.
! 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;
!                   Kod >=2: se tiparesc si (marimile/pas): y,z,w si f(y),f(z),f(w)si numitorii la z, w.
!       : -IESIRE : cod de terminare/eroare:
!		 0:	toleranta "eps" atinsa;
!		 1:	depasirea nr.maxim de iteratii admis "lnit";
!        2: |f(x)| <=eps_f
!        2: Proces stationar;
!		-1: divergenta
!       -2:	Numitor < eps_num
!       -3,-4,-5: y = x; z = y; w = z;
!
logical print_pas
logical stationar
integer width
logical IsInfinity

real num    ! numitorul in formule

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

x =x0

!eps_num =1E-12

! Pas 0: x
call compute_f(x, fx)
! Test fx: for "Infinity" (overflow)
IF(IsInfinity(fx)) THEN
write(1,*) '* fx = Infinit'
kod =-10
RETURN
ENDIF
! Test fx: for "eps_f"
if(abs(fx) <=eps_f) then
kod =2
goto 500
endif

call compute_f1(x, f1x)
!*
if(print_pas) then
write(1,*); write(1,*) 'iter:', iter;
write(1,*) 'x, fx:', x, fx
write(1,*) 'f1x:', f1x
endif

IF(IsInfinity(f1x)) THEN
write(1,*) '* f1x = Infinit'
kod =-10
RETURN
ENDIF

num =f1x
in =check_num(num, 'f''(x)')
    if(in <0) then
    kod =-2
    goto 500
    endif

! Pas 1: y
y =x -fx/f1x

IF(isnan(y)) THEN
write(1,*) '* y = NaN'
kod =-10
RETURN
ENDIF

call compute_f(y, fy)
if(print_pas) then
write(1,*) 'y, fy:', y,fy
endif

! Test fy:
IF(IsInfinity(fy)) THEN
write(1,*) '* fy = Infinit'
kod =-10
RETURN
ENDIF

if(abs(fy) <=eps_f) then
kod =2; x =y;
goto 500
endif

!Test y-x:
if(y ==x) then
kod =-3
    goto 500
endif

! Test numitor pentru z:
num =(fx -2*fy)
in =check_num(num, 'Numitor_z')
    if(in <0) then
    kod =-2
    goto 500
    endif

z =y -fx/(fx -2*fy)*fy/f1x

IF(isnan(z)) THEN
write(1,*) '* z = NaN'
kod =-10
RETURN
ENDIF

call compute_f(z, fz)
if(print_pas) then
write(1,*) 'z, fz:', z,fz
write(1,*) 'num-z:', num
endif

IF(IsInfinity(fz)) THEN
write(1,*) '* fz = Infinit'
kod =-10
RETURN
ENDIF

if(abs(fz) <=eps_f) then
kod =2; x =z;
goto 500
endif

if(z ==y) then
kod =-4
    goto 500
endif

num =div_dif(x,z)*div_dif(y,z)
in =check_num(num, 'Numitor_w')
    if(in <0) then
    kod =-2
    goto 500
    endif

w =z -fz*div_dif(x,y)/div_dif(x,z)/div_dif(y,z)*(1+ fz/fx)
IF(isnan(w)) THEN
write(1,*) '* w = NaN'
kod =-10
RETURN
ENDIF

call compute_f(w, fw)
!**
if(print_pas) then
write(1,*) 'w, fw:', w,fw
write(1,*) 'num-w:', num
endif
!*
IF(IsInfinity(fw)) THEN
write(1,*) '* fw = Infinit'
kod =-10
RETURN
ENDIF
if(abs(fw) <=eps_f) then
kod =2; x =w;
goto 500
endif


if(w ==z) then
kod =-5
    goto 500
endif

!* Check dif_div_2:
d21 =div_dif_2(y,x,z)
IF(isnan(d21)) THEN
write(1,*) '* f[y,x,z] = NaN'
kod =-10
RETURN
ENDIF

d22 =div_dif_2(y,x,w)
IF(isnan(d22)) THEN
write(1,*) '* f[y,x,w] = NaN'
kod =-10
RETURN
ENDIF

d23 =div_dif_2(z,x,w)
IF(isnan(d23)) THEN
write(1,*) '* f[z,x,w] = NaN'
kod =-10
RETURN
ENDIF

num =div_dif(x,w) +( div_dif_2(y,x,z) -div_dif_2(y,x,w) -div_dif_2(z,x,w) )*(x-w)
in =check_num(num, 'Numitor_x')
    if(in <0) then
    kod =-2
    goto 500
    endif

dif_w =fw/num

!x =w -fw/( div_dif(x,w) +( div_dif_2(y,x,z) -div_dif_2(y,x,w) -div_dif_2(z,x,w) )*(x-w) )
x =w -dif_w

if(print_pas) then
write(1,*) 'num-x:', num
write(1,*) 'x, dif, x:', x, dif_w
endif

IF(IsInfinity(x)) THEN
write(1,*) '* x = Infinit'
kod =-10
!RETURN
ENDIF

IF(isnan(x)) THEN
write(1,*) '* x = NaN'
kod =-10
!RETURN
ENDIF

500 continue
END
