  We discuss the renormalisation group (RG) evolution for the $\Delta S = 1$
operators in unquenched QCD with $N_f = 3$ ($m_u=m_d=m_s$) or, more generally,
$N_f = 2+1$ ($m_u=m_d \ne m_s$) flavors. In particular, we focus on the
specific problem of how to treat the singularities which show up only for
$N_f=3$ or $N_f = 2+1$ in the original solution of Buras {\it et al.} for the
RG evolution matrix at next-to-leading order. On top of Buras {\it et al.}'s
original treatment, we use a new method of analytic continuation to obtain the
correct solution in this case. It is free of singularities and can therefore be
used in numerical analysis of data sets calculated in lattice QCD.
