docs: 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