  We present results from two different methods for the calculation of hadron
spectra in QCD and SUSY QCD with large primary energies $\sqrt{s}$ up to
$10^{16}$ GeV. The two methods considered are a Monte Carlo (MC) simulation and
the evolution of fragmentation functions described by the
Dokshitzer-Gribov-Lipatov-Altarelli-Parisi (DGLAP) equations. We find that the
pion, nucleon and all-hadron spectra calculated with the two methods agree
well. The MC simulation is performed with new hadronization functions (in
comparison with our previous work), motivated by low-energy ($\sqrt s < M_Z$)
data and DGLAP. The hadron spectra calculated with both sets of hadronization
functions agree well, which indicates that our method for calculating the
hadronization function works successfully. The small difference in the
calculated hadron spectra characterizes the uncertainties of this method. We
calculate also the spectra of photons, neutrinos and nucleons and compare them
with other published results. The agreement is good for all $x$ from $\sim
10^{-5}$ up to $x\leq 0.3$. The consistency of the spectra calculated by
different methods allows to consider the spectral shape as a signature of
models with decays or annihilations of superheavy particles, such as
topological defects or superheavy DM. The UHECR spectra from these sources are
calculated.
