The first step is the completion of the whole constraint list (primary, secondary , ternary etc.), and checking that no secondary constraint leads to a contradiction (i.e., empty constraint surface).
Remark: The time evolution of the constraints is performed with the "total" Hamiltonian:
$H_T(p,q) = p\dot{q}-L+\sum_{\alpha}\lambda_{\alpha} \phi^{(1)}_{\alpha} $
where $\phi^{(1)}_{\alpha}$ are the primary constraints and $\lambda_{\alpha}$ are Lagrange multipliers.
The next step is the computation of the Poisson bracket matrix of all constraints:
$P_{\alpha\beta} = \{\phi_{\alpha}, \phi_{\beta}\}|_{\Sigma}$
(for all constraints $\phi_1, ., ., ., \phi_n$). $\Sigma$ is the constraint surface.
The number of first class constraints is equal to the corank of the matrix $P$: $n-\mathrm{rank}(P)$.