Molecular dynamics computer simulation methods are very important for understanding mechanisms of chemical, physical, and biological processes. The reliability of molecular dynamics simulations strongly depends on the integration schemes used in the simulations. In this work, we developed new rigid body integration schemes for molecular dynamics simulations. Our approach is based on a numerically exact solution to the free rigid body problem, which is used in the classical propagator splitting scheme. We use the Taylor series expansion of rotational dynamical variables in conjunction with the recursive solution for higher order derivatives of these variables. Such an approach is computationally very efficient, robust, and easy to implement, and it does not employ Jacobi elliptic functions, while still providing the numerically exact solution of the free rigid body problem. Our studies showed that the new integration methods have long-time stability and accuracy properties which are comparable to those of existing symplectic integrators. The extension to the case of a canonical ensemble is also developed, allowing one to perform simulations at constant temperatures.