The operator canonical perturbation theory (CPT) is an efficient tool for solving the molecular vibration-rotation Schrödinger equation. The corresponding Watson Hamiltonian can be written using angular momentum cylindrical ladder operators (Jz, J± = Jx ∓ iJy) possessing the Lie algebra su(2) commutation relations [J+, J-] = 2Jz, [Jz, J±] = ±J±. The reduced effective Hamiltonians suitable for fitting to observed spectra are traditionally based on Hermitian basis sets, e.g., (J2)lJz n,(J+ m+J- m)+. It is beneficial to re-express such Hamiltonians using sums of normal ordered products of ladder operators Jz aJ+ bJ- c. For instance, in the CPT, the unitary transformations reduce to the normal ordering problem. Similarly, the dipole moment transition probabilities can be evaluated using Wigner functions, D0,ε 1(ε=0,±1), possessing complex commutation relations with Jα-operators. The related line strengths are proportional to the squared matrix elements of the unitary transformed dipole moment operator given by a polynomial in normal ordered products D0,ε 1Jz aJ+ bJ- c. We have applied the technique of combinatorial calculations associated with the classical representation theory and obtained compact formulas for normal ordering of the typical raw products Jz aJ+ bJ- cJz dJ+ eJ- f and Jz aJ+ bJ- cD0,ε 1. The theory of universal enveloping Lie algebras and representation theory show that the resulting formulas cannot be further improved from the mathematical point of view. The results have a wide scope of applications in other fields of physics, including other operators with the same commutation properties. The resulting expressions for normal orderings and routines for evaluation of matrix elements, as well as the systematic tests of working formulas, were coded in Fortran and the software is available through the GitHub hosting service.