  Dimensional regularization is incompatible with the standard covariant
projection methods that are used to calculate the short-distance coefficients
in inclusive heavy quarkonium production and annihilation rates. A new method
is developed that allows dimensional regularization to be used consistently to
regularize the infrared and ultraviolet divergences that arise in these
perturbative calculations. We illustrate the method by calculating the leading
color-octet terms and the leading color-singlet terms in the gluon
fragmentation functions for arbitrary quarkonium states. We resolve a
discrepancy between two previous calculations of the gluon fragmentation
functions for the spin-triplet P-wave quarkonium states.
