def _equilibrium_payoffs_abreu_sannikov(rpg, tol=1e-12, max_iter=500,
                                        u_init=np.zeros(2)):
    """
    Using 'abreu_sannikov' algorithm to compute the set of payoff pairs
    of all pure-strategy subgame-perfect equilibria with public randomization
    for any repeated two-player games with perfect monitoring and
    discounting, following Abreu and Sannikov (2014).

    Parameters
    ----------
    rpg : RepeatedGame
        Two player repeated game.

    tol : scalar(float), optional(default=1e-12)
        Tolerance for convergence checking.

    max_iter : scalar(int), optional(default=500)
        Maximum number of iterations.

    u_init : ndarray(float, ndim=1), optional(default=np.zeros(2))
        The initial guess of threat points.

    Returns
    -------
    hull : scipy.spatial.ConvexHull
        The convex hull of equilibrium payoff pairs.

    References
    ----------
    .. [1] Abreu, Dilip, and Yuliy Sannikov. "An algorithm for
       two‐player repeated games with perfect monitoring." Theoretical
       Economics 9.2 (2014): 313-338.
    """
    sg, delta = rpg.sg, rpg.delta

    if sg.N != 2:
        msg = "this algorithm only applies to repeated two-player games."
        raise NotImplementedError(msg)

    best_dev_gains = _best_dev_gains(rpg)
    IC = np.empty(2)
    action_profile_payoff = np.empty(2)
    # auxiliary array for checking if payoff is inside the convex hull
    # first two entries for payoff point, and the last entry is 1.
    extended_payoff = np.ones(3)
    # array to store new points of C in each intersection
    # at most 4 new points will be generated
    new_pts = np.empty((4, 2))
    # array to store the points of W
    # the length of v is limited by |A1|*|A2|*4
    W_new = np.empty((np.prod(sg.nums_actions)*4, 2))
    W_old = np.empty((np.prod(sg.nums_actions)*4, 2))
    # count the new points generated in each iteration
    n_new_pt = 0

    # copy the threat points
    u = np.copy(u_init)

    # initialization
    payoff_pts = \
        sg.payoff_profile_array.reshape(np.prod(sg.nums_actions), 2)
    W_new[:np.prod(sg.nums_actions)] = payoff_pts
    n_new_pt = np.prod(sg.nums_actions)

    n_iter = 0
    while True:
        W_old[:n_new_pt] = W_new[:n_new_pt]
        n_old_pt = n_new_pt
        hull = ConvexHull(W_old[:n_old_pt])

        W_new, n_new_pt = \
            _R(delta, sg.nums_actions, sg.payoff_arrays,
               best_dev_gains, hull.points, hull.vertices,
               hull.equations, u, IC, action_profile_payoff,
               extended_payoff, new_pts, W_new)

        n_iter += 1
        if n_iter >= max_iter:
            break

        # check convergence
        if n_new_pt == n_old_pt:
            if np.linalg.norm(W_new[:n_new_pt] - W_old[:n_new_pt]) < tol:
                break

        # update threat points
        _update_u(u, W_new[:n_new_pt])

    hull = ConvexHull(W_new[:n_new_pt])

    return hull