module defpi
    implicit none
    real(16),parameter::pi=4.0Q0*atan(1.0Q0)
end module defpi

program f90_dka
!***********************************************************
! DKA method (Durand-Kerner-Aberth)
! f(x)=a0*x^n-a1*x^(n-1)+....+an=0
! http://www.ecs.shimane-u.ac.jp/~kyoshida/c6(2002).pdf
!***********************************************************
    use defpi
    implicit none
    integer,parameter::itmax=10000
    real(16),parameter::eps=1.0Q-20
    real(16),parameter::zero=1.0Q-1000
    integer::i,j,n,iii
    real(16),allocatable::a(:),z1(:),z2(:),err(:)
    real(16),allocatable::ff0(:),z1mem(:),z2mem(:)
    integer,allocatable::imem(:),icoef(:)
    real(16)::w0,w1,w2,f1,f2,ffe
    real(16)::p1,p2,ee,ae,a1,a2,sume
    real(16)::rr,rc,theta,xc,yc
    character(len=50)::fnameR
    character(len=100)::str

!----------------------------------------------------
!   Data input
!----------------------------------------------------
    call getarg(1,fnameR)
    open(11,file=fnameR,status='old')
    read(11,'(a)') str
    read(11,*) n
    allocate(a(0:n),z1(1:n),z2(1:n),err(1:n))
    allocate(ff0(1:n),z1mem(1:n),z2mem(1:n))
    allocate(imem(1:n),icoef(1:n))
    do i=0,n
        read(11,*) a(i)
    end do
    close(11)
!----------------------------------------------------
!   Treatment of coefficient
!----------------------------------------------------
    w0=a(0)
    do i=0,n
        a(i)=a(i)/w0
    end do
!----------------------------------------------------
!   Initial value of x
!----------------------------------------------------
    write(6,'(a)') trim(adjustl(str))
!    rc=RINI1(n,a)
    rc=RINI2(n,a)
    do i=1,n
        theta=2.0Q0*pi/dble(n)*(dble(i)-0.75Q0)
        z1(i)=-a(1)/dble(n)+rc*cos(theta)
        z2(i)=rc*sin(theta)
    end do
!----------------------------------------------------
!   Output of initial value
!----------------------------------------------------
    ! Write to the standard output
    write(6,'(a)') 'Initial value'
    write(6,'(a)') '    rc             -a(1)/n'
    write(6,'(2(e15.7))') rc,-a(1)/dble(n)
    write(6,'(a)') '    i    z1             z2             rr'
    do i=1,n
        rr=sqrt((z1(i)+a(1)/dble(n))**2+z2(i)**2)
        write(6,'(i5,3(e15.7))') i,z1(i),z2(i),rr
    end do
    ! Write to the file for GMT
    open(12, file='_iteration.txt', status='replace')
    do i=1,n
        write(12,'(2(e15.7))') z1(i),z2(i)
    end do
!----------------------------------------------------
!   Iteration by DK method
!----------------------------------------------------
    do i=1,n
        ff0(i)=1.0Q1000
    end do
    do iii=1,itmax
        ee=0.0Q0
        sume=0.0Q0
        do i=1,n
            ! Numerator
            f1=1.0Q0
            f2=0.0Q0
            do j=1,n
                w1=f1*z1(i)-f2*z2(i)
                w2=f2*z1(i)+f1*z2(i)
                f1=w1+a(j)
                f2=w2
            end do
            ! Denominator
            p1=1.0Q0
            p2=0.0Q0
            do j=1,n
                if(j==i)cycle
                w1=p1*(z1(i)-z1(j))-p2*(z2(i)-z2(j))
                w2=p1*(z2(i)-z2(j))+p2*(z1(i)-z1(j))
                p1=w1
                p2=w2
            end do
            a1=(f1*p1+f2*p2)/(p1*p1+p2*p2)
            a2=(f2*p1-f1*p2)/(p1*p1+p2*p2)
            ! Evaluation of (f/df)
            ae=sqrt(a1*a1+a2*a2)
            err(i)=ae
            if(ee<ae)ee=ae
            z1(i)=z1(i)-a1
            z2(i)=z2(i)-a2
            if(abs(z1(i))<zero)z1=0.0Q0
            if(abs(z2(i))<zero)z2=0.0Q0
            ! Evaluation of (f)
            f1=1.0Q0
            f2=0.0Q0
            do j=1,n
                w1=f1*z1(i)-f2*z2(i)
                w2=f2*z1(i)+f1*z2(i)
                f1=w1+a(j)
                f2=w2
            end do
            ffe=sqrt(f1*f1+f2*f2)
            if(ffe<ff0(i))then
                ff0(i)=ffe
                imem(i)=iii
                z1mem(i)=z1(i)
                z2mem(i)=z2(i)
            end if
        end do
        do i=1,n
            write(12,'(2(e15.7))') z1(i),z2(i)
        end do
        sume=sume+ae
        if(ee<eps)exit
        if(iii==itmax)exit
    end do
    close(12)
!----------------------------------------------------
!   Output of result
!----------------------------------------------------
    do i=1,n
        icoef(i)=i
    end do
    call BISORT(n,z1,icoef)
    ! Write to the standard output
    write(6,'(a)') '    iii Average error'
    write(6,'(i7,e15.7)') iii,sume/dble(n)
    write(6,'(a)') '    i    z1             z2             rr             error'
    do i=1,n
        j=icoef(i)
        rr=sqrt((z1(j)+a(1)/dble(n))**2+z2(j)**2)
        write(6,'(i5,4(e15.7))') j,z1(j),z2(j),rr,err(j)
    end do
    write(6,'(a)') 'Values at minimum f(x)'
    write(6,'(a)') '    i    z1             z2             f(z)           iii'
    do i=1,n
        j=icoef(i)
        write(6,'(i5,3(e15.7),i7)') j,z1mem(j),z2mem(j),ff0(j),imem(j)
    end do
    ! Write to the file for GMT
    call GMTDAT(n,iii,rc,-a(1)/dble(n),0.0Q0)

stop

contains
!----------------------------------------------------
!   Sort of result
!----------------------------------------------------
   subroutine BISORT(n,z1,icoef)
        ! Sorting in small order */
        integer,intent(in)::n
        real(16),intent(in)::z1(1:n)
        integer,intent(inout)::icoef(1:n)
        integer::i,j,iwork
        real(16)::wz(1:n),work

        do i=1,n
            wz(i)=z1(i)
        end do
        do i=1,n
            do j=i+1,n
                if(wz(j)<wz(i))then
                    work=wz(i)
                    wz(i)=wz(j)
                    wz(j)=work

                    iwork=icoef(i)
                    icoef(i)=icoef(j)
                    icoef(j)=iwork
                end if
            end do
        end do
    end subroutine BISORT
!----------------------------------------------------
!   Data output for GMT drawing
!----------------------------------------------------
    subroutine GMTDAT(n,iii,rr,zxc,zyc)
        integer,intent(in)::n,iii
        real(16),intent(in)::rr,zxc,zyc
        integer::i,k
        real(16)::zz1(0:iii,1:n),zz2(0:iii,1:n)
        real(16)::zx,zy

        open(11,file='_iteration.txt',status='old')
        do k=0,iii
            do i=1,n
                read(11,*) zz1(k,i),zz2(k,i)
            end do
        end do
        close(11)

        open(12, file='_gmt_inp.txt', status='replace')
        do i=1,n
            write(12,'(a)') trim(adjustl('>'))
            do k=0,iii
                write(12,'(2(e15.7))') zz1(k,i),zz2(k,i)
            end do
        end do
        close(12)

        open(12, file='_data_circ.txt', status='replace')
        do i=0,360
            zx=zxc+rr*cos(dble(i)/180.0Q0*pi)
            zy=zyc+rr*sin(dble(i)/180.0Q0*pi)
            write(12,'(2(e15.7))') zx,zy
        end do
        close(12)

    end subroutine GMTDAT
!----------------------------------------------------
!   Calculation of r-value for initial x-value
!   (simple method)
!----------------------------------------------------
    real(16) function RINI1(n,a)
        integer,intent(in)::n
        real(16),intent(in)::a(0:n)
        integer::i,ic
        real(16)::rm

        ic=0
        do i=1,n
            if(eps<abs(a(i)))ic=ic+1
        end do
        rm=0.0D0
        do i=1,n
            rr=dble(ic)*abs(a(i))**(1.0Q0/dble(i))
            if(rm<rr)rm=rr
        end do
        RINI1=rm
    end function RINI1
!----------------------------------------------------
!   Calculation of r-value for initial x-value
!   (General method)
!----------------------------------------------------
    real(16) function RINI2(n,a)
        integer,intent(in)::n
        real(16),intent(in)::a(0:n)

        integer,parameter::iimax=100
        real(16),parameter::eps=1.0Q-6
        integer::i,ii
        real(16)::b(0:n)
        real(16)::xx,dx,ff,df

        ! Set of coefficient
        do i=0,n
            b(i)=-abs(a(i))
        end do
        b(0)=abs(b(0))
        b(1)=0.0Q0
        ! Search initial x-value
        xx=0.0Q0
        dx=1.0Q0
        do
            xx=xx+dx
            ff=b(0)
            do i=1,n
                ff=xx*ff+b(i)
            end do
            if(0.0Q0<ff)exit
        end do
        ! Iteration of Newton method
        write(6,'(a)') 'Radius for initial value'
        write(6,'(a)') '   ii    x              f(x)'
        do ii=1,iimax
            ff=b(0)
            do i=1,n
                ff=xx*ff+b(i)
            end do
            df=n*b(0)
            do i=1,n-1
                df=xx*df+(n-i)*b(i)
            end do
            xx=xx-ff/df
            write(6,'(i5,2(e15.7))') ii,xx,ff
            if(abs(ff/df)<eps)exit
        end do
        ! Answer
        RINI2=xx
    end function RINI2
!----------------------------------------------------
!   End of program
!----------------------------------------------------
end program f90_dka
