  We revisit the challenging problem of finding an efficient Monte Carlo (MC)
algorithm solving the constrained evolution equations for the initial-state QCD
radiation. The type of the parton (quark, gluon) and the energy fraction x of
the parton exiting emission chain (entering hard process) are predefined, i.e.
constrained throughout the evolution. Such a constraint is mandatory for any
realistic MC for the initial state QCD parton shower. We add one important
condition: the MC algorithm must not require the a priori knowledge of the full
numerical exact solutions of the evolution equations, as is the case in the
popular ``Markovian MC for backward evolution''. Our aim is to find at least
one solution of this problem that would function in practice. Finding such a
solution seems to be definitely within the reach of the currently available
computer CPUs and the sophistication of the modern MC techniques. We describe
in this work the first example of an efficient solution of this kind. Its
numerical implementation is still restricted to the pure gluon-strahlung. As
expected, it is not in the class of the so-called Markovian MCs. For this
reason we refer to it as belonging to a class of non-Markovian MCs. We show
that numerical results of our new MC algorithm agree very well (to 0.2%) with
the results of the other MC program of our own (unconstrained Markovian) and
another non-MC program QCDnum16. This provides a proof of the existence of the
new class of MC techniques, to be exploited in the precision perturbative QCD
calculations for the Large Hadron Collider.
