Once, we numerically integrated the precession and nutation of a spheroidal rigid Earth (Kubo and Fukushima 1987). As a natural extension, we tried to integrate the rotation of a triaxial rigid Earth numerically and faced a problem: a loss of precision in long-term integration. This is due to the smallness of the characteristic period of the problem: 1 day. Of course, one can integrate the rotational motion in higher precision arithmetics with a smaller stepsize. However, the quadruple precision integration is roughly 30 times more time-consuming than the double precision integration. See Table 1. Therefore, it is desirable if there is a formulation 1) reducing the overall integration error, 2) being independent on the choice of the integrator and 3) requiring no extra computations. The key points to achieve this goal will be to find a set of variables which 1) are efficiently convertible to the physical quantities required finally, say, the orientation matrix in the case of the rotational dynamics, and 2) vary with time as smoothly as possible. In this note, we report a discovery of such an example.