anchor.diff:  
     function sort_complex_conjugate_pairs(roots) BIND(C)
       implicit none
  -    integer :: sort_complex_conjugate_pairs
  -    double complex, dimension(4), intent(inout) :: roots
  -    integer :: npairs = 0
  -    integer :: i, j
  -    double complex :: r, r_conj
  +    integer(c_int) :: sort_complex_conjugate_pairs
  +    complex(c_double_complex), dimension(4), intent(inout) :: roots
  +    integer(c_int) :: npairs, i, j
  +    complex(c_double_complex) :: r, r_conj
  +    npairs = 0
       i = 1
       do while (i <= 3)
          r = roots(i)
          if (abs(DIMAG(r)) > IMAG_TOLERANCE) then
  -          do j = 1, 3-i
  +          do j = 1, 4-i
                r_conj = roots(i+j)
                if (       abs(DREAL(r) - DREAL(r_conj)) < ZERO_TOLERANCE &
                     .and. abs(DIMAG(r) + DIMAG(r_conj)) < ZERO_TOLERANCE) then
  
positive.code:        endif
         S = 0.5*sqrt(SSx4)
      endif
      sqrtp = sqrt(-SSx4 - 2*p + q/3)
      sqrtm = sqrt(-SSx4 - 2*p - q/3)
      out(1) = -0.25*b - S + 0.5*sqrtp;
      out(2) = -0.25*b - S - 0.5*sqrtp;
      out(3) = -0.25*b + S + 0.5*sqrtm;
  
negative.code: MODULE poly_solvers
    USE iso_c_binding, only: c_double, c_double_complex, c_int
    IMPLICIT NONE
  
    real(c_double), parameter :: ZERO_TOLERANCE = 1.d-15
    real(c_double), parameter :: IMAG_TOLERANCE = 1.d-14
    real(c_double), parameter :: IMAG_TOLERANCE_SQRD = IMAG_TOLERANCE**2
    real(c_double), parameter :: PI = acos(-1.d0)
    real(c_double), parameter :: PIx2 = 2*PI
    complex(c_double_complex), dimension(3), parameter :: CUBE_ROOTS_OF_1 = (/ cmplx(1.d0,0.d0, kind=c_double),        &
                                                                               exp(cmplx(0.d0,PIx2/3, kind=c_double)), &
                                                                               exp(cmplx(0.d0,2*PIx2/3, kind=c_double)) /)
  
  CONTAINS
  
    SUBROUTINE quartic_solve (Poly, out) BIND(C)
      implicit none
      real(c_double), dimension(5), intent(in) :: Poly
      complex(c_double_complex), dimension(4), intent(out) :: out
      real(c_double) :: e, d, c, b, bb, p, q, r, cc, dd, ee, ccc, Delta, Dee, Delta_0, Delta_1
      complex(c_double_complex) :: S, SSx4, phi, zQ, SSx4_max, sqrtp, sqrtm
      real(c_double) :: abs_SSx4_max, abs_SSx4
      integer(c_int) :: ir
      e = Poly(1) / Poly(5)
      d = Poly(2) / Poly(5)
      c = Poly(3) / Poly(5)
      b = Poly(4) / Poly(5)
      bb = b*b
      p = 0.125d0 * (8*c - 3*bb)
      q = 0.125d0 * (bb*b - 4*b*c + 8*d)
      r = (-3*bb*bb + 256*e - 64*b*d + 16*bb*c) / 256
      c = p
      d = q
  