  We present the first calculation at next-to-leading order (NLO) in $\alpha_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\epsilon$ dimensions. We also extract the divergences in the virtual NLO
corrections analytically, and detail the cancellation of all divergences after
renormalization. The NLO corrections have a dramatic effect on the shape of the
fragmentation function, and they significantly increase the fragmentation
probability.
