A QCD factorization formalism was recently proposed for evaluating heavy quarkonium production at large pT at collider energies. With systematically calculated short-distance partonic hard parts and evolution kernels of fragmentation functions (FFs), the predictive power of this factorization approach relies on our knowledge of a large number of universal FFs at an input factorization scale μ0≳2mQ with heavy quark mass mQ. With the large heavy quark mass, the relative motion of the heavy quark pair inside a heavy quarkonium is effectively nonrelativistic. We evaluate these universal input FFs using nonrelativistic QCD (NRQCD) factorization and express the large number of FFs in terms of a few universal NRQCD long-distance matrix elements with perturbatively calculated coefficients. We derive complete contributions to the single-parton FFs at both O(αs) and O(αs2) and the heavy quark-pair FFs at O(αs). We present detailed derivation for all contributions involving long-distance matrix elements of S-wave NRQCD QQ¯ states (P-wave contributions in a companion paper, Y.-Q. Ma, J.-W. Qiu, and H. Zhang, arXiv:1401.0524 [Phys. Rev. D (to be published)]). Our results bridge the gap between the QCD factorization formalism and its phenomenological applications.