We present the first calculation at next-to-leading order (NLO) in α s of a fragmentation function into quarkonium whose form at leading order is a nontrivial function of z , namely the fragmentation function for a gluon into a spin-singlet S-wave state at leading order in the relative velocity. To calculate the real NLO corrections, we introduce a new subtraction scheme that allows the phase-space integrals to be evaluated in 4 dimensions. We extract all ultraviolet and infrared divergences in the real NLO corrections analytically by calculating the phase-space integrals of the subtraction terms in 4 − 2 ϵ dimensions. We also extract the divergences in the virtual NLO corrections analytically, and detail the cancellation of all divergences...