  Exploiting recently obtained analytic solutions of classical Yang-Mills
equations for higher order perturbations in the field of the dilute object
(proton), we derive the complete first saturation correction to the single
inclusive semi-hard gluon production in high energy proton-nucleus collisions
by applying the Lehmann-Symanzik-Zimmermann reduction formula. We thus finalize
the program started by Balitsky (see Ref.~\cite{Balitsky:2004rr}) and
independently by Chirilli, Kovchegov and Wertepny (see
Ref.~\cite{Chirilli:2015tea}) albeit using a very different approach to carry
out our calculations. We extracted the functional dependence of gluon spectrum
on the color charge densities of the colliding objects; thus our results can be
used to evaluate complete first saturation correction to the double/multiple
inclusive gluon productions.
