Based on a rigorous relativistic treatment of rotational motion within the first post-Newtonian approximation of Einstein's theory of gravity, we have developed a new numerical theory of rigid Earth rotation. This theory is constructed by numerical integration and treats for the first time the whole spectrum of relativistic issues in a consistent way: (1) relativistic time scales, (2) relativistic scaling of astronomical constants, (3) relativistic torques and (4) geodetic precession as an additional torque in the equations of rotational motion. We will discuss the fundamental ideas of the theory and compare our results with the currently used models. While in the quasi-Newtonian limit our theory reproduces SMART97 within the accuracy of the latter, we have found that the standard procedure to treat geodetic precession leads to errors in all modern theories of Earth rotation with a magnitude of up to 200 ?as within 100 years. The possibility to extend our numerical code to model the rotational motion of other solar system bodies will be addressed as well.