anchor.diff:      next_candidates(num_next_candidates)%width = half_width
   
     end subroutine split_candidate
  +
  +  subroutine allocate_candidates(num_candidates, next_candidates)
  +    integer(c_int), intent(in) :: num_candidates
  +    type(LocateCandidate), allocatable, intent(inout) :: next_candidates(:)
  +
  +    if (allocated(next_candidates)) then
  +       if (size(next_candidates) < 4 * num_candidates) then
  +          deallocate(next_candidates)
  +          allocate(next_candidates(4 * num_candidates))
  +       end if
  +    else
  +       allocate(next_candidates(4 * num_candidates))
  +    end if
  +
  +  end subroutine allocate_candidates
  +
  +  subroutine update_candidates( &
  +       num_nodes, degree, point, num_candidates, candidates, &
  +       num_next_candidates, next_candidates)
  +
  +    integer(c_int), intent(in) :: num_nodes, degree
  +    real(c_double), intent(in) :: point(1, 2)
  +    integer(c_int), intent(in) :: num_candidates
  +    type(LocateCandidate), intent(in) :: candidates(:)
  +    integer(c_int), intent(inout) :: num_next_candidates
  +    type(LocateCandidate), allocatable, intent(inout) :: next_candidates(:)
  +    ! Variables outside of signature.
  +    integer(c_int) :: cand_index
  +    logical(c_bool) :: predicate
  +
  +    call allocate_candidates(num_candidates, next_candidates)
  +
  +    do cand_index = 1, num_candidates
  +       call contains_nd( &
  +            num_nodes, 2, candidates(cand_index)%nodes, &
  +            point(1, :), predicate)
  +
  +       if (predicate) then
  +          num_next_candidates = num_next_candidates + 4
  +          call split_candidate( &
  +               num_nodes, degree, candidates(cand_index), &
  +               num_next_candidates, next_candidates)
  +       end if
  +    end do
  +
  +  end subroutine update_candidates
   
     subroutine locate_point( &
          num_nodes, nodes, degree, x_val, y_val, s_val, t_val) &
  
positive.code:     ! We divide by `3n` rather than `n` since we have tracked thrice the
      ! centroid rather than the centroid itself (to avoid round-off until
      ! right now).
      s_approx = ( &
           sum(candidates(:num_candidates)%centroid_x) / &
           (3.0_dp * num_candidates))
      t_approx = ( &
           sum(candidates(:num_candidates)%centroid_y) / &
           (3.0_dp * num_candidates))
  
      call newton_refine( &
           num_nodes, nodes, degree, x_val, y_val, &
  
negative.code:     real(c_double), intent(out) :: s_val, t_val
      ! Variables outside of signature.
      real(c_double) :: point(1, 2)
      integer(c_int) :: sub_index
      integer(c_int) :: num_candidates, num_next_candidates
      type(LocateCandidate), allocatable :: candidates_odd(:), candidates_even(:)
      real(c_double) :: s_approx, t_approx
      real(c_double) :: actual(1, 2)
      logical(c_bool) :: is_even
  
      point(1, :) = [x_val, y_val]
      ! Start out with the full curve.
      allocate(candidates_odd(1))  ! First iteration is odd.
      candidates_odd(1) = LocateCandidate(1.0_dp, 1.0_dp, 1.0_dp, nodes)
      num_candidates = 1
      s_val = LOCATE_MISS
  
      is_even = .TRUE.  ! At zero.
      do sub_index = 1, MAX_LOCATE_SUBDIVISIONS + 1
         is_even = .NOT. is_even  ! Switch parity.
         num_next_candidates = 0
  
         if (is_even) then
            ! Read from even, write to odd.
            call update_candidates( &
                 num_nodes, degree, point, num_candidates, candidates_even, &
                 num_next_candidates, candidates_odd)
         else
            ! Read from odd, write to even.
            call update_candidates( &
                 num_nodes, degree, point, num_candidates, candidates_odd, &
                 num_next_candidates, candidates_even)
         end if
  
         ! If there are no more candidates, we are done.
         if (num_next_candidates == 0) then
            return
         end if
  
         num_candidates = num_next_candidates
  
      end do
  