! Simpla precizie
! Metoda Laguerre pentru polinoame
! Zero-uri complexe
! n: gradul polinomului.
! coef: vectorul coef. incepand cu gradul cel mai mare
! rad: tablou complex pt. radacini
!
! Polinom: p(z) =a_n*x^n +a_(n-1)*x^(n-1) + ... +a_1*z +a_0
! a_i = reali
!
! Fisier de date:
! 1) Titlul problemei: <=80 caractere.
! 2) Grad polinom; numar de radacini; cod tiparire corectii; cod limba rezultate: n, n_rad, kod_print, lang
!   character(1) lang: r, R =romana; alt char. =engleza; kod_print =0/#0: Nu/Da: tiparste iteratiile.
! 3) Coeficientii (reali): a_n, a_(n-1), ..., a_0
! 4) eps, nlim: toleranta; nr, limita de iteratii.
! 5) x0(1:n_rad): aproximatiile initiale
! 6) [Optional: lang: text, primul caracter: r,R/alt char. =romana/engleza]

use comun

implicit real(8) (a-h, o-z)
character(260) fname, dummy*80
!real(4), allocatable :: coef(:), coef_1(:), coef_2(:)
!coef(0:n); coef_1(0:n-1); coef_2(0:n-2)

!complex, allocatable :: rad(:)
!complex, allocatable :: x0(:)
!integer, allocatable :: iter(:)

complex(8) z01, rad1
complex(8) zpolval, polval, polval_i
logical is_Complex
integer status, width
character(1) lang
character(100) z_char
character(50) end_text(0:2,2), pol_val(2), adj_eps(2)

! Ending - initialization:
end_text(0,1) ='Toleranta atinsa.'
end_text(1,1) ='* LNIT depasit.'
end_text(2,1) ='* Proces stationar.'
!end_text(-1)=' * Divergenta.'

end_text(0,2) ='Tolerance met.'
end_text(1,2) ='* LNIT exceeded'
end_text(2,2) ='* Stationary Process.'
!end_text(-1)=' * Divergence.'

pol_val(1) =' Valori polinom'
pol_val(2) =' Polynomial Values'

adj_eps(1) =': Ajustare EPS:'
adj_eps(2) =': Ajusting EPS:'


call GetFile(fname, status, 'Select Input File'C)
open(1, file =fname, status ='old', form ='formatted')

read(1,'(a)') dummy		! Title
read (1, *) n, n_rad, kod_print

if(n_rad >n .or. n_rad ==0) n_rad =n

allocate (coef(0:n))
read (1,*) (coef(i), i=n,0,-1)

print *, 'Coeficienti: ', (coef(i), i=n,0,-1)

allocate (rad(1:n_rad), z0(1:n_rad))
allocate (iter(1:n_rad), kod(1:n_rad))
iter =0

read (1,*) eps, nlim
read (1,*) (z0(1:n_rad))

lang =' '
read (1,*, err=110, end =110) lang

rom = lang =='r' .or. lang =='R'
en =.not.rom
goto 120

110 en =.true.

120 continue
if(rom) then
ilang =1
else
ilang =2
endif

! Coef. derivatelor:
allocate (coef_1(0:n-1), coef_2(0:n-2))

do i = n-1, 0, -1
coef_1(i) =coef(i+1)*(i+1)
enddo

do i = n-2, 0, -1
coef_2(i) =coef_1(i+1)*(i+1)
enddo

if(kod_print /=0) write(1,*) 

! Apel LAG:
do i_rad =1, n_rad
iter1 =nlim
z01 =z0(i_rad)

! Adjust EPS:
eps_min =ulp(abs(z01))
eps_i =eps

if(eps_i < eps_min) then
eps_i =eps_min
write(1, '(/, I2, a, 1pE9.2)') i_rad, trim(adj_eps(ilang)), eps_min
endif


call LAG (z01, eps_i, iter1, rad1, kod1, kod_print)
iter(i_rad) =iter1
rad(i_rad) =rad1
kod(i_rad) =kod1
enddo

! Roots:
if(rom) then
write(1,'(/,a, t27,a, t43,a)') ' Radacini', 'Iteratii', 'Terminare'
else
write(1,'(/,a, t27,a, t43,a)') ' Roots', 'Iterations', 'Ending'
endif

do i_rad =1,n_rad
rad1 =rad(i_rad)
iter1 =iter(i_rad); kod1 =kod(i_rad);
iwiter =width(iter1)

if(is_Complex(rad1)) then
write(z_char,*) rad1
write(1, '(I2, a, 4x,I<iwiter>, t43,a)') i_rad, trim(z_char), iter(i_rad), end_text(kod1, ilang)

else
write(z_char,*) real(rad1)
write(1, '(I2, a, t27,I<iwiter>, t43,a)') i_rad, trim(z_char), iter(i_rad), end_text(kod1, ilang)
endif
enddo

! Polynomial values:
write(1,'(/, a)') pol_val(ilang)

do i_rad =1,n_rad
polval_i =zpolval(n, coef, rad(i_rad))
rad1 =rad(i_rad)
if(is_Complex(rad1)) then
write(z_char,*) polval_i
write(1, '(I2, a)') i_rad, trim(z_char)

else
write(z_char,*) real(polval_i)
write(1, '(I2, a)') i_rad, trim(z_char)
endif

enddo


END

