Numerical integrations in celestial mechanics often involve the repeated computation of a rotation with a constant angle. A direct evaluation of these rotations yields a linear drift of the distance to the origin. This is due to roundoff in the representation of the sine $s$ and cosine $c$ of the angle $\theta$. In a computer, one generally gets $c^2 + s^2 \ne 1$, resulting in a mapping that is slightly contracting or expanding. In the present paper we present a method to find pairs of representable real numbers $s$ and $c$ such that $c^2 + s^2$ is as close to 1 as possible. We show that this results in a drastic decrease of the systematic error, making it negligible compared to the random error of other operations. We also verify that this approach gives good results in a realistic celestial mechanics integration.