roots Subroutine

   pure subroutine roots(a, na, croots, irl)
      real(dp), intent(in) :: a(:) !! Array of polynomial coefficients (lowest power first)
      integer(ip), intent(in) :: na !! Number of coefficients
      real(dp), intent(inout) :: croots(:,:) !! Initial guesses and output roots (real, imaginary)
      integer(ip), intent(in) :: irl !! Number of roots to solve for

      integer(ip) :: i, ir, n1, n2, j, iter
      real(dp) :: z(2), zs(2), cb(2), cc(2), cb_old(2), dif, denom, temp
      real(dp), parameter :: conv_tol = 1.0e-14_dp
      real(dp), parameter :: abs_tol = 1.0e-15_dp
      integer(ip), parameter :: max_iter = 100

      ir = 0
      n1 = na - 1
      n2 = n1 - 1

      do while (ir < irl)
         z(1) = croots(ir+1,1)
         z(2) = croots(ir+1,2)

         iter = 0
         do
            iter = iter + 1
            if (iter > max_iter) then
               ! Exceeded maximum iterations - exit with current approximation
               exit
            end if

            cb(1) = a(n1+1)
            cb(2) = 0.0_dp
            cc(1) = a(n1+1)
            cc(2) = 0.0_dp

            do i = 0, n2
               j = n2 - i
               ! Save old cb values before updating
               cb_old(1) = cb(1)
               cb_old(2) = cb(2)
               ! Update cb (polynomial value)
               temp = (z(1) * cb(1) - z(2) * cb(2)) + a(j+1)
               cb(2) = z(1) * cb(2) + z(2) * cb(1)
               cb(1) = temp
               ! Update cc (derivative) using OLD cb values
               if (j /= 0) then
                  temp = (z(1) * cc(1) - z(2) * cc(2)) + cb_old(1)
                  cc(2) = (z(1) * cc(2) + z(2) * cc(1)) + cb_old(2)
                  cc(1) = temp
               end if
            end do

            zs(1) = z(1)
            zs(2) = z(2)

            ! Newton's method
            denom = cc(1) * cc(1) + cc(2) * cc(2)
            if (denom < 1.0e-30_dp) then
               ! Derivative is zero - cannot continue Newton's method
               exit
            end if
            z(1) = z(1) - ((cb(1) * cc(1) + cb(2) * cc(2)) / denom)
            z(2) = z(2) + ((cb(1) * cc(2) - cb(2) * cc(1)) / denom)

            ! Convergence criterion with protection against division by zero
            if (abs(zs(1)) > abs_tol) then
               dif = abs((zs(1) - z(1)) / zs(1))
            else
               ! Use absolute difference when zs(1) is near zero
               dif = abs(zs(1) - z(1))
            end if

            if (abs(zs(2)) > abs_tol) then
               dif = dif + abs((zs(2) - z(2)) / zs(2))
            else
               ! Use absolute difference when zs(2) is near zero
               dif = dif + abs(zs(2) - z(2))
            end if

            if (dif <= conv_tol) exit
         end do

         croots(ir+1,1) = z(1)
         croots(ir+1,2) = z(2)
         ir = ir + 1
      end do

   end subroutine roots