To compute the class of states for the given transition probability matrix I have been given the following transition probability matrix of a markov chain:
$P = \begin{pmatrix}
\frac{3}{4}{} & 0 & \frac{1}{4} &0 \\ 
\frac{1}{2} & 0 & 0 & \frac{1}{2}\\ 
 0&  0&  0& 1\\ 
 \frac{1}{4}&\frac{1}{4}  &\frac{1}{4}  &\frac{1}{4} 
\end{pmatrix}$
I want to compute the class of each state. My current understanding of class of state is this:
The class of state $i$, denoted by, $C(i)$ is the set of all such states which are communicating with $i$. Where communicating means, two states will be communicating if they can be accessed from each other, irrespective of the number of steps, meaning that one $i$ state can be accessed through another state $j$ in $m$ steps and $j$ can be accessed through $i$ in $n$ number of steps. $m$ and $n$ need not be equal.
Please let me know if my current understanding of class is the correct one. Then, to actually compute the class, I draw the digraph. I don't know how to do that in latex. So, I am uploading an image. Here, we go:

From this digraph, I want to compute the class of states $1,2,3,$ and $4$.
$C(1) = (2,4)$
1 is communicating with 2 although we cannot go from 1 to 2 directly but we can go like this $1 => 3 => 4 => 2$. Also, we can come to 2 from 1 like this: $1=>3=>4=>2$. Is this the correct way to solve these type of problems?
 A: This is correct.  It is, however, nice to have a numerical algorithm to solve the problem when the number of states grows large.  There's a simple one.
The matrix $\mathbb A = (a_{ij})$ having ones where $\mathbb P$ is nonzero and otherwise having zeros is the adjacency matrix of the Markov chain.  That is, for any states $i$ and $j,$ $a_{ij}=1$ means there is a nonzero chance of a  transition from $i$ to $j.$
When you multiply $\mathbb A$ by itself using the usual rules of matrix multiplication but interpret $+$ as Boolean "or" and $*$ as Boolean "and," the result is the adjacency matrix $\mathbb{A}^2$ describing all two-step transitions.  That is, $(\mathbb{A}^2)_{ij}=1$ if and only if there is positive probability of a transition from $i$ to $j$ that takes two steps.  Consequently,
$$\mathbb{A}^{[2]} = \mathbb{A} + \mathbb{A}^{[1]}\,\mathbb{A}$$
(where $\mathbb{A}^{[1]} = \mathbb{A}$) describes all transitions that can occur within at most two steps.  Inductively,
$$\mathbb{A}^{[n+1]} = \mathbb{A} + \mathbb{A}^{[n]}\,\mathbb{A}$$
is the adjacency matrix for all transitions that can occur within at most $n+1$ steps.
Eventually the sequence of $\mathbb{A}^{[n]}$ stabilizes because although a $0$ entry can turn into a $1,$ a $1$ entry can never revert to $0.$  This occurs when $n$ is the length of the longest distance between any two states, which will happen by the time $n$ equals the number of states.  Thus, this sequence of powers of $\mathbb A$ has a limit $$\mathbb{A}^{*} = \lim_{n\to\infty} \mathbb{A}^{[n]}$$ which can be directly computed in at most $n-1$ steps.
The communicating classes can be read off $\mathbb{A}^{*}:$ a state $i$ communicates with a state $j$ if and only if $(\mathbb{A}^{*})_{ij}=1.$

To implement this algorithm, you can use the ordinary arithmetic of positive integers.  When $x,y\in\{0,1\},$ define
$$x+y = \min(x + y, 1)$$
and otherwise let $xy$ be computed as usual.  Here is an R implementation as a function that takes the transition matrix $\mathbb P$ as input (assumed to be square, of course) and returns $\mathbb{A}^{*}.$
star <- function(P) {
  A.star <- A <- ifelse(P != 0, 1, 0)
  for (i in seq_len(dim(A)[1])) {
    A.new <- pmin(A.star + A.star %*% A, 1)
    if(isTRUE(all.equal(A.new, A.star))) break
    A.star <- A.new
  }
  A.star
}

For example, with P <- matrix(c(3,2,0,1,  0,0,0,1,  1,0,0,1,  0,2,4,1), 4) / 4 as in the question, the function call
star(P)

produces the output

     [,1] [,2] [,3] [,4]
[1,]    1    1    1    1
[2,]    1    1    1    1
[3,]    1    1    1    1
[4,]    1    1    1    1


indicating all classes communicate.
As an example where there are separate communicating classes, consider a Markov chain on five states where $1$ stays fixed, $2$ and $3$ transition to each other with probability $1/2,$ and $4$ and $5$ transition to each other with probability $1/2.$  Obviously they comprise three communicating classes $\{1\},$ $\{2,3\},$ and $\{4,5\}.$  Here is the code and its output:
P <- matrix(c(1,0,0,0,0, 0,0,1/2,0,0, 0,1/2,0,0,0, 0,0,0,0,1/2, 0,0,0,1/2,0), 5)
(A <- star(P))


     [,1] [,2] [,3] [,4] [,5]
[1,]    1    0    0    0    0
[2,]    0    1    1    0    0
[3,]    0    1    1    0    0
[4,]    0    0    0    1    1
[5,]    0    0    0    1    1


The classes are straightforward to find, as this R implementation shows.  Its output is a list whose entries are arrays of indexes into $\mathbb{A}$ (starting at $1$) to represent the states.
classes <- function(A) {
  k <- seq_len(dim(A)[1])
  C <- list()
  classid <- 0
  while(length(k) > 0) {
    C[[classid <- classid+1]] <- j <- which(A[k[1],]==1)
    k <- setdiff(k, j)
  }
  C
}

For example, the call
classes(A)

(using the preceding output) produces

[[1]]
[1] 1

[[2]]
[1] 2 3

[[3]]
[1] 4 5


listing the three classes as arrays of indexes into A.
