def _all_intersections(nodes_first, nodes_second):
    r"""Find the points of intersection among a pair of curves.

    .. note::

       There is also a Fortran implementation of this function, which
       will be used if it can be built.

    .. note::

       This assumes both curves are in :math:`\mathbf{R}^2`, but does not
       **explicitly** check this. However, functions used here will fail if
       that assumption fails, e.g. :func:`bbox_intersect` and
       :func:`newton_refine() <._intersection_helpers._newton_refine>`.

    Args:
        nodes_first (numpy.ndarray): Control points of a curve to be
            intersected with ``nodes_second``.
        nodes_second (numpy.ndarray): Control points of a curve to be
            intersected with ``nodes_first``.

    Returns:
        Tuple[numpy.ndarray, bool]: An array and a flag:

        * A ``2 x N`` array of intersection parameters.
          Each row contains a pair of values :math:`s` and :math:`t`
          (each in :math:`\left[0, 1\right]`) such that the curves
          intersect: :math:`B_1(s) = B_2(t)`.
        * Flag indicating if the curves are coincident.

    Raises:
        ValueError: If the subdivision iteration does not terminate
            before exhausting the maximum number of subdivisions.
        NotImplementedError: If the subdivision process picks up too
            many candidate pairs. This typically indicates tangent
            curves or coincident curves (though there are mitigations for
            those cases in place).
    """
    curve_first = SubdividedCurve(nodes_first, nodes_first)
    curve_second = SubdividedCurve(nodes_second, nodes_second)
    candidate1 = Linearization.from_shape(curve_first)
    candidate2 = Linearization.from_shape(curve_second)
    # Handle the line-line intersection case as a one-off.
    both_linear, result = check_lines(candidate1, candidate2)
    if both_linear:
        return result

    candidates = [(candidate1, candidate2)]
    intersections = []
    coincident = False
    for _ in six.moves.xrange(_MAX_INTERSECT_SUBDIVISIONS):
        candidates = intersect_one_round(candidates, intersections)
        if len(candidates) > _MAX_CANDIDATES:
            candidates = prune_candidates(candidates)
            # If pruning didn't fix anything, we check if the curves are
            # coincident and "fail" if they aren't.
            if len(candidates) > _MAX_CANDIDATES:
                params = coincident_parameters(nodes_first, nodes_second)
                if params is None:
                    raise NotImplementedError(
                        _TOO_MANY_TEMPLATE.format(len(candidates))
                    )

                else:
                    intersections = params
                    coincident = True
                    # Artificially empty out candidates so that this
                    # function exits.
                    candidates = []
        # If none of the candidate pairs have been accepted, then there are
        # no more intersections to find.
        if not candidates:
            if intersections:
                # NOTE: The transpose of a C-ordered array is Fortran-ordered,
                #       i.e. this is on purpose.
                return np.array(intersections, order="C").T, coincident

            else:
                return np.empty((2, 0), order="F"), coincident

    msg = _NO_CONVERGE_TEMPLATE.format(_MAX_INTERSECT_SUBDIVISIONS)
    raise ValueError(msg)