This work presents an elegant formalism to model the evolution of the full two rigid body problem. The equations of motion, given in a Cartesian coordinate system, are expressed in terms of spherical harmonics and Wigner D-matrices. The algorithm benefits from the numerous recurrence relations satisfied by these functions allowing a fast evaluation of the mutual potential. Moreover, forces and torques are straightforwardly obtained by application of ladder operators taken from the angular momentum theory and commonly used in quantum mechanics. A numerical implementation of this algorithm is made. Tests show that the present code is significantly faster than those currently available in literature.