  We investigate a large class of perturbative QCD (pQCD) renormalization
schemes whose beta functions $\beta(a)$ are meromorphic functions of the
running coupling and give finite positive value of the coupling $a(Q^2)$ in the
infrared regime ("freezing"), $a(Q^2) \to a_0$ for $Q^2 \to 0$. Such couplings
automatically have no singularities on the positive axis of the squared momenta
$Q^2$ ($ \equiv -q^2$). Explicit integration of the renormalization group
equation (RGE) leads to the implicit (inverted) solution for the coupling, of
the form $\ln (Q^2/Q^2_{\rm in}) = {\cal H}(a)$. An analysis of this solution
leads us to an algebraic algorithm for the search of the Landau singularities
of $a(Q^2)$ on the first Riemann sheet of the complex $Q^2$-plane, i.e., poles
and branching points (with cuts) outside the negative semiaxis. We present
specific representative examples of the use of such algorithm, and compare the
found Landau singularities with those seen after the 2-dimensional numerical
integration of the RGE in the entire first Riemann sheet, where the latter
approach is numerically demanding and may not always be precise. The specific
examples suggest that the presented algebraic approach is useful to find out
whether the running pQCD coupling has Landau singularities and, if yes, where
precisely these singularities are.
