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: MODULE poly_solvers
    USE iso_c_binding, only: c_double
    IMPLICIT NONE
  
    real(c_double), bind(C, name="R") :: R = 0.02625
    double precision, parameter :: ZERO_TOLERANCE = 1.d-15
    double precision, parameter :: IMAG_TOLERANCE = 1.d-14
    double precision, parameter :: IMAG_TOLERANCE_SQRD = IMAG_TOLERANCE**2
    double precision, parameter :: PI = acos(-1.d0)
    double precision, parameter :: PIx2 = 2*PI
    double precision, dimension(3), parameter :: z_loc = (/ 0.d0, 1.d0, 0.d0 /)
    double complex, dimension(3), parameter :: CUBE_ROOTS_OF_1 = (/ cmplx(1.d0,0.d0),        &
                                                                    exp(cmplx(0.d0,PIx2/3)), &
                                                                    exp(cmplx(0.d0,2*PIx2/3)) /)
  
  CONTAINS
  
    SUBROUTINE quartic_solve (Poly, out) BIND(C)
      implicit none
      double precision, dimension(5), intent(in) :: Poly
      double complex, dimension(4), intent(out) :: out
      double precision :: e, d, c, b, bb, p, q, r, cc, dd, ee, ccc, Delta, Dee, Delta_0, Delta_1
      double complex :: S, SSx4, phi, zQ, SSx4_max, sqrtp, sqrtm
      double precision :: abs_SSx4_max, abs_SSx4
      integer :: 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.125 * (8*c - 3*bb)
      q = 0.125 * (bb*b - 4*b*c + 8*d)
      r = (-3*bb*bb + 256*e - 64*b*d + 16*bb*c) / 256
      c = p
      d = q
  
negative.code:        endif
         S = 0.5*sqrt(SSx4)
      endif
      sqrtp = sqrt(-SSx4 - 2*p + q/S)
      sqrtm = sqrt(-SSx4 - 2*p - q/S)
      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;
  