  Within the Quark-Gluon String Model A.B. Kaidalov found a behaviour of quark
and diquark fragmentation functions for $z\rightarrow 0$ and $z\rightarrow 1$,
and proposed interpolation formulae for the functions in the whole region of
$z$. These functions must be a solution of the well-known system of the
integral equations. A simplified Monte Carlo estimation of the functions, based
on usage of the fragmentation functions at $z\rightarrow 1$ as the kernel
functions of the system, does not reproduce Kaidalov's results. An improvement
of the Monte Carlo simulations is proposed in this paper. It can be implemented
in Monte Carlo event generators such as Los Alamos QGSM, QGSJet-II and the
Geant4 QGS model. It will improve a description of experimental data in the
models, especially, the description of the latest NA61/SHINE Collaboration data
on $\pi{\rm C}$ interactions. Description of the data is a problem in DPMJet,
QGSJet, EPOS and Sibyll models.
