! Simpla precizie
! POLINOM PERTURBAT
! p~(x) = p(x) +eps*q(x)
!
! Metoda Laguerre pentru polinoame
! + K = nr. cond. pentru stabilitatea radacinilor
!
! 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*x +a_0
! a_i = reali
!
! Polinom: q(z) =b_(nq)*x^n +b_(nq-1)*x^(nq-1) + ... +b_1*x +b_0
! b_i = reali
!
! Fisier de date:
! 1) Titlul problemei: <=80 caractere.
! 2) Grad polinom "p"; cod tiparire corectii: n, kod_print
! 3) Coeficientii lui "p" (reali): a_n, a_(n-1), ..., a_0
! 4) toleranta; nr. limita de iteratii: tol, nlim;
! 5) aproximatiile initiale - rad. lui "p": x00(1:n_rad)
! perturbare:
! 6) Grad "q(x)"; perturbare: n_q, eps
! 7) Coeficientii "q": b(n_q), ..., b(0)
!
! kod_print: 0/#0: tiparire corectii
!             >=2: tiparire K
!
! Daca: n_q =0 sau eps =0:
! - Se calculeaza numai radacinile lui "p"
! - Linia (7) de date poate lipsi.
!
use work

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 :: x00(:), x0(:)
!complex, allocatable :: rad(:)
!integer, allocatable :: iter(:)

complex x01, rad1, dif
complex polval, polval_i, z
integer status
complex K1, zK1

logical is_real


call GetFile(fname, status, 'Selecteaza Fisierul de Intrare'C)

open(1, file =fname, status ='old', form ='formatted')

read(1,'(a)') dummy		! Title

read (1,*) n, kod_print	! gradul

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

allocate (x00(1:n))

read (1,*) tol, nlim
read (1,*) x00(1:n)

! "q":
read (1,*) n_q, eps

IF(n_q ==0 .or. eps ==0.) GOTO 100

allocate (coef_q(0:n))  ! acelasi cu "n"
coef_q =0.

read (1,*) (coef_q(i), i=n_q,0,-1)


100 print *, 'Coeficienti (P): '
print *, (coef_0(i), i=n,0,-1)

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

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

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

allocate (x0(1:n))
allocate (iter(1:n))
iter =0

! I) Rad. lui "P":
! ---------------
write(1,'(/, a)') 'Radacinile lui "P"'
write(1,'(a)')    '------------------'

conjig =.false.
rad1 =0.

P_DO: DO i_rad =1, n
iter1 =nlim
x01 =x00(i_rad)

! Adjust tol:
tol_min =ulp(abs(x01))
tol_i =tol

if(tol_i < tol_min) then
tol_i =tol_min
write(1, '(/,I2, a, 1pE9.2)') i_rad, ': Ajustare tol:', tol_min
endif

call LAG (x01, tol_i, iter1, i_rad, rad1, kod_print)
150 continue
iter(i_rad) =iter1
x0(i_rad) =rad1

ENDDO P_DO

! Roots (P):
write(1,'(/, a, t43,a)') ' Radacini:', 'Iteratii'
do i_rad =1,n
rad1 =x0(i_rad)

if(is_real(rad1)) then
    write(1,*) i_rad, real(rad1), iter(i_rad)
else
    write(1,*) i_rad, rad1, iter(i_rad)
endif

enddo

! II) Rad. lui "P~":
! -----------------
IF(n_q ==0 .or. eps ==0.) THEN
! Polynomial "P" values:
write(1,'(/, a)') ' Valori polinom:'
do i_rad =1,n
polval_i =polval(n, coef, x0(i_rad))
    if(is_real(polval_i)) then
write(1,*), i_rad, real(polval_i)
    else
write(1,*), i_rad, polval_i
    endif
enddo

GOTO 500
ENDIF

! Coef. lui "p-perturbat":
coef =coef +EPS*coef_q

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

! Coef. derivatelor - Initializare:
coef_1(0:n-1) =0.
coef_2(0:n-2) =0.

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

allocate (rad(1:n))
rad =0.
iter =0

write(1,'(/, a)') 'Radacinile lui "P-perturbat"'
write(1,'(a)')    '---------------------------'

! Apel LAG:
rad1 =0.

Pq_DO: DO i_rad =1, n
iter1 =nlim
x01 =x0(i_rad)

! Adjust tol:
tol_min =ulp(abs(x01))
tol_i =tol

if(tol_i < tol_min) then
tol_i =tol_min
write(1, '(/,I2, a, 1pE9.2)') i_rad, ': Ajustare tol:', tol_min
endif

call LAG (x01, tol_i, iter1, i_rad, rad1, kod_print)
250 continue
iter(i_rad) =iter1
rad(i_rad) =rad1

ENDDO Pq_DO

! Roots (P~):
write(1,'(/, a, t43,a)') ' Radacini:', 'Iteratii'
do i_rad =1,n
rad1 =rad(i_rad)

if(is_real(rad1)) then
    write(1,*) i_rad, real(rad1), iter(i_rad)
    else
    write(1,*) i_rad, rad1, iter(i_rad)
    endif

enddo

! Polynomial values:
write(1,'(/, a)') ' Valori polinom:'
do i_rad =1,n
polval_i =polval(n, coef, rad(i_rad))
    if(is_real(polval_i)) then
write(1,*), i_rad, real(polval_i)
    else
write(1,*), i_rad, polval_i
    endif
enddo

! K:
Print_K: IF(kod_print >=2) THEN
write(1,'(/,a, t30,a, t47,a, t62,a)') '  Radacina', 'K1(z)', 'K1(z)*eps', 'z(eps)-z'     !'z -K1*eps'

do i_rad =1,n
rad1 =rad(i_rad)
dif =rad1 -x0(i_rad)
zK1 =K1(rad1)

if(is_real(rad1)) then
    write(1, '(1pG14.7, t28,1pG14.6, t45,1pG11.3, t60,1pG11.3)') real(rad1), real(zK1), real(zK1*eps), real(dif)
else
    write(1, '( (1pG14.7,1pG14.7), t28,1pG14.6,a1, t45,1pG11.3,a1, t60,1pG11.3,a1,: 2x,a)') &
    rad1, abs(zK1),'*', abs(zK1*eps), '*', abs(dif), '*'
endif

enddo

write(1,'(/,a)') ' * in modul'
write(1,'(a)') ' ( K1(z) =-q(z)/p''(z) )'
ENDIF Print_K

500 continue
END
