Implement the imitation game algorithm by McLennan and Tourky (2006)
    for computing an approximate fixed point of `T`.

    Parameters
    ----------
    is_approx_fp : callable
        A callable with signature `is_approx_fp(v)` which determines
        whether `v` is an approximate fixed point with a bool return
        value (i.e., True or False)

    For the other parameters, see Parameters in compute_fixed_point.

    Returns
    -------
    x_new : scalar(float) or ndarray(float)
        Approximate fixed point.

    converged : bool
        Whether the routine has converged.

    iterate : scalar(int)
        Number of iterations.

def _compute_fixed_point_ig(T, v, max_iter, verbose, print_skip, is_approx_fp,
                            *args, **kwargs):
    """
    Implement the imitation game algorithm by McLennan and Tourky (2006)
    for computing an approximate fixed point of `T`.

    Parameters
    ----------
    is_approx_fp : callable
        A callable with signature `is_approx_fp(v)` which determines
        whether `v` is an approximate fixed point with a bool return
        value (i.e., True or False)

    For the other parameters, see Parameters in compute_fixed_point.

    Returns
    -------
    x_new : scalar(float) or ndarray(float)
        Approximate fixed point.

    converged : bool
        Whether the routine has converged.

    iterate : scalar(int)
        Number of iterations.

    """
    if verbose == 2:
        start_time = time.time()
        _print_after_skip(print_skip, it=None)

    x_new = v
    y_new = T(x_new, *args, **kwargs)
    iterate = 1
    converged = is_approx_fp(x_new)

    if converged or iterate >= max_iter:
        if verbose == 2:
            error = np.max(np.abs(y_new - x_new))
            etime = time.time() - start_time
            print_skip = 1
            _print_after_skip(print_skip, iterate, error, etime)
        if verbose >= 1:
            if not converged:
                warnings.warn(_non_convergence_msg, RuntimeWarning)
            elif verbose == 2:
                print(_convergence_msg.format(iterate=iterate))
        return x_new, converged, iterate

    if verbose == 2:
        error = np.max(np.abs(y_new - x_new))
        etime = time.time() - start_time
        _print_after_skip(print_skip, iterate, error, etime)

    # Length of the arrays to store the computed sequences of x and y.
    # If exceeded, reset to min(max_iter, buff_size*2).
    buff_size = 2**8
    buff_size = min(max_iter, buff_size)

    shape = (buff_size,) + np.asarray(x_new).shape
    X, Y = np.empty(shape), np.empty(shape)
    X[0], Y[0] = x_new, y_new
    x_new = Y[0]

    tableaux = tuple(np.empty((buff_size, buff_size*2+1)) for i in range(2))
    bases = tuple(np.empty(buff_size, dtype=int) for i in range(2))
    max_piv = 10**6  # Max number of pivoting steps in _lemke_howson_tbl

    while True:
        y_new = T(x_new, *args, **kwargs)
        iterate += 1
        converged = is_approx_fp(x_new)

        if converged or iterate >= max_iter:
            break

        if verbose == 2:
            error = np.max(np.abs(y_new - x_new))
            etime = time.time() - start_time
            _print_after_skip(print_skip, iterate, error, etime)

        try:
            X[iterate-1] = x_new
            Y[iterate-1] = y_new
        except IndexError:
            buff_size = min(max_iter, buff_size*2)
            shape = (buff_size,) + X.shape[1:]
            X_tmp, Y_tmp = X, Y
            X, Y = np.empty(shape), np.empty(shape)
            X[:X_tmp.shape[0]], Y[:Y_tmp.shape[0]] = X_tmp, Y_tmp
            X[iterate-1], Y[iterate-1] = x_new, y_new

            tableaux = tuple(np.empty((buff_size, buff_size*2+1))
                             for i in range(2))
            bases = tuple(np.empty(buff_size, dtype=int) for i in range(2))

        m = iterate
        tableaux_curr = tuple(tableau[:m, :2*m+1] for tableau in tableaux)
        bases_curr = tuple(basis[:m] for basis in bases)
        _initialize_tableaux_ig(X[:m], Y[:m], tableaux_curr, bases_curr)
        converged, num_iter = _lemke_howson_tbl(
            tableaux_curr, bases_curr, init_pivot=m-1, max_iter=max_piv
        )
        _, rho = _get_mixed_actions(tableaux_curr, bases_curr)

        if Y.ndim <= 2:
            x_new = rho.dot(Y[:m])
        else:
            shape_Y = Y.shape
            Y_2d = Y.reshape(shape_Y[0], np.prod(shape_Y[1:]))
            x_new = rho.dot(Y_2d[:m]).reshape(shape_Y[1:])

    if verbose == 2:
        error = np.max(np.abs(y_new - x_new))
        etime = time.time() - start_time
        print_skip = 1
        _print_after_skip(print_skip, iterate, error, etime)
    if verbose >= 1:
        if not converged:
            warnings.warn(_non_convergence_msg, RuntimeWarning)
        elif verbose == 2:
            print(_convergence_msg.format(iterate=iterate))

    return x_new, converged, iterate