  In the framework of the nonrelativistic quantum chromodynamics factorization
formalism, we study the processes of $\psi(nS)$ and $\Upsilon(nS)$ decay into a
lepton pair or a charm pair associated with two jets up to the next-to-leading
order in velocity expansion. We present the analytic expressions for the
differential decay rate to the invariant mass of the lepton pair or charm pair.
We find that the ratio of the next-to-leading order short-distance coefficient
to the leading order one is in the range from -5.5 to -12.4. The relativistic
corrections are so large that they modify the leading order prediction
significantly. Utilizing the analytic expressions, we also investigate the
relativistic corrections in different kinematic regions and their dependence on
the masses of the initial-state quarkonium and the final-state fermion. In
addition, we study the momentum distribution of $D^{*+}$ in the process
$\Upsilon(1S)\to c\bar{c}gg\to D^{*+}X$.
