We study the splitting functions for the evolution of fragmentation distributions and the coefficient functions for single-hadron production in semi-inclusive e + e − annihilation in massless perturbative QCD for small values of the momentum fraction and scaling variable x, where their fixed-order approximations are completely destabilized by huge double logarithms of the form α n s x −1 ln 2n−a x. Complete analytic all-order expressions in Mellin-N space are presented for the resummation of these terms at the next-to-next-to-leading logarithmic accuracy. The poles for the first moments, related to the evolution of hadron multiplicities, and the small-x instabilities of the next-to-leading order splitting and coefficient functions are removed by this resummation, which leads to an oscillatory small-x behaviour and functions that can be used at N = 1 and down to extremely small values of x. First steps are presented towards extending these results to the higher accuracy required for an all-x combination with the state-of-the-art next-to-next-to-leading order large-x results.If the constants up to F (m) n,ℓ are known for all n and ℓ, then the splitting functions and coefficient functions can be determined at N m LL accuracy at all orders of the strong coupling. As observed in Ref. [19], the n th order small-x contributions to F T, φ are built up from n terms of the form(2.4)Since the terms with ε −2n+1 , . . . , ε −n−1 have to cancel in sum (2.1), there are n−1 relations between the LL coefficients A n,k which lead to the constants F (0) n,ℓ in Eq. (2.3), n−2 relations between the NLL coefficients B n,k etc. As discussed above, a N m LO calculation fixes the (nonvanishing) coefficients of ε −n , . . . , ε −n+m at all orders n, adding m + 1 more relations between the coefficients in Eq. (2.4). Consequently the highest m+1 double logarithms, i.e., the N m LL approximation, can be determined order by order from the N m LO results. Finally the resulting series, here calculated to order α 18 s using FORM and TFORM [24], can be employed to find their generating functions via over-constrained systems of linear equations. The whole procedure is analogous to, if computationally more involved than, the large-x resummation in Ref. [25].