Find all the intersection points between the current convex hull
    and the two IC constraints. It is done by iterating simplex
    counterclockwise.

    Parameters
    ----------
    C : ndarray(float, ndim=2)
        The 4 by 2 array for storing the generated potential
        extreme points of one action profile. One action profile
        can only generate at most 4 points.

    points : ndarray(float, ndim=2)
        Coordinates of the points in the W, which construct a
        feasible payoff convex hull.

    vertices : ndarray(float, ndim=1)
        Indices of points forming the vertices of the convex hull,
        which are in counterclockwise order.

    equations : ndarray(float, ndim=2)
        [normal, offset] forming the hyperplane equation of the facet

    extended_payoff : ndarray(float, ndim=1)
        The array [payoff0, payoff1, 1] for checking if
        [payoff0, payoff1] is in the feasible payoff convex hull.

    IC : ndarray(float, ndim=1)
        The minimum IC continuation values.

    tol : scalar(float)
        The tolerance for checking if two values are equal.

    Returns
    -------
    C : ndarray(float, ndim=2)
        The generated potential extreme points.

    n : scalar(int)
        The number of found intersection points.

def _find_C(C, points, vertices, equations, extended_payoff, IC, tol):
    """
    Find all the intersection points between the current convex hull
    and the two IC constraints. It is done by iterating simplex
    counterclockwise.

    Parameters
    ----------
    C : ndarray(float, ndim=2)
        The 4 by 2 array for storing the generated potential
        extreme points of one action profile. One action profile
        can only generate at most 4 points.

    points : ndarray(float, ndim=2)
        Coordinates of the points in the W, which construct a
        feasible payoff convex hull.

    vertices : ndarray(float, ndim=1)
        Indices of points forming the vertices of the convex hull,
        which are in counterclockwise order.

    equations : ndarray(float, ndim=2)
        [normal, offset] forming the hyperplane equation of the facet

    extended_payoff : ndarray(float, ndim=1)
        The array [payoff0, payoff1, 1] for checking if
        [payoff0, payoff1] is in the feasible payoff convex hull.

    IC : ndarray(float, ndim=1)
        The minimum IC continuation values.

    tol : scalar(float)
        The tolerance for checking if two values are equal.

    Returns
    -------
    C : ndarray(float, ndim=2)
        The generated potential extreme points.

    n : scalar(int)
        The number of found intersection points.
    """
    n = 0
    weights = np.empty(2)
    # vertices is ordered counterclockwise
    for i in range(len(vertices)-1):
        n = _intersect(C, n, weights, IC,
                       points[vertices[i]],
                       points[vertices[i+1]], tol)

    n = _intersect(C, n, weights, IC,
                   points[vertices[-1]],
                   points[vertices[0]], tol)

    # check the case that IC is an interior point of the convex hull
    extended_payoff[:2] = IC
    if (np.dot(equations, extended_payoff) <= tol).all():
        C[n, :] = IC
        n += 1

    return C, n