docs: def time_relaxations_direct(P, p0, obs, times=[1]):
      r"""Compute time-relaxations of obs with respect of given initial distribution.
  
      relaxation(k) = p0 P^k obs
  
      Parameters
      ----------
      P : ndarray, shape=(n, n) or scipy.sparse matrix
          Transition matrix
      p0 : ndarray, shape=(n)
          initial distribution
      obs : ndarray, shape=(n)
          Vector representing observable on discrete states.
      times : array-like, shape(n_t)
          Vector of time points at which the (auto)correlation will be evaluated
  
      Returns
      -------
      relaxations : ndarray, shape(n_t)
      """
      n_t = len(times)
      times = np.sort(times)
  
      # maximum time > number of rows?
      if times[-1] > P.shape[0]:
          use_diagonalization = True
          R, D, L = rdl_decomposition(P)
          # discard imaginary part, if all elements i=0
          if not np.any(np.iscomplex(R)):
              R = np.real(R)
          if not np.any(np.iscomplex(D)):
              D = np.real(D)
          if not np.any(np.iscomplex(L)):
              L = np.real(L)
          rdl = (R, D, L)
  
      f = np.empty(n_t, dtype=D.dtype)
  
      if use_diagonalization:
          for i in range(n_t):
              f[i] = time_relaxation_direct_by_diagonalization(P, p0, obs, times[i], rdl)
      else:
          start_values = None
          for i in range(n_t):
              f[i], start_values = time_relaxation_direct_by_mtx_vec_prod(P, p0, obs, times[i], start_values, True)
      return f