  The standard analytic solution of the renormalization group (RG) evolution
for the $\Delta S = 1$ Wilson coefficients involves several singularities,
which complicate analytic solutions. In this paper we derive a singularity-free
solution of the next-to-leading order (NLO) RG equations, which greatly
facilitates the calculation of $\epsilon_K^{\prime}$, the measure of direct
$CP$ violation in $K\to \pi\pi$ decays. Using our new RG evolution and the
latest lattice results for the hadronic matrix elements, we calculate the ratio
$\epsilon_{K}^{\prime}/\epsilon_{K}$ (with $\epsilon_{K}$ quantifying indirect
$CP$ violation) in the Standard Model (SM) at NLO to
$\epsilon_{K}^{\prime}/\epsilon_{K} = (1.06 \pm 5.07) \times 10^{-4} $, which
is $2.8\,\sigma$ below the experimental value. We also present the evolution
matrix in the high-energy regime for calculations of new physics contributions
and derive easy-to-use approximate formulae. We find that the RG amplification
of new-physics contributions to Wilson coefficients of the electroweak penguin
operators is further enhanced by the NLO corrections: If the new contribution
is generated at the scale of 1-10 TeV, the RG evolution between the new-physics
scale and the electroweak scale enhances these coefficients by 50-100 %. Our
solution contains a term of order $\alpha_{EM}^2/\alpha_s^2$, which is
numerically unimportant for the SM case but should be included in studies of
high-scale new-physics.
