In the framework of the nonrelativistic quantum chromodynamics factorization formalism, we study the processes of psi (nS) and gamma (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 gamma(1S) -> c (c) over bar gg -> D*+X.