We compare two approaches to evaluate cross sections of heavy-quarkonium production at next-to-leading order in nonrelativistic QCD involving $S$- and $P$-wave Fock states: the customary approach based on phase space slicing and the approach based on dipole subtraction recently elaborated by us. We find reasonable agreement between the numerical results of the two implementations, but the dipole subtraction implementation outperforms the phase space slicing one both with regard to accuracy and speed.