Motivated by detecting Earth twins in extreme precise radial velocity data, PEXO is a package aiming for extremely high-precision modeling of the motion of single stars or stars in a system. PEXO is general enough to account for binary motion and stellar reflex motions induced by dark companions and is precise enough to treat various relativistic effects both in the solar system and in the target system. We also model the post-Newtonian barycentric motion for future tests of general relativity in binaries. We benchmark PEXO with the pulsar timing package TEMPO2 and find that PEXO produces numerically similar results with timing precision of about 1 nanosecond, space-based astrometry to a precision of 1 μas, and radial velocity of 1 μm/s and improves on TEMPO2 for decade-long timing data of nearby targets, due to its consideration of third-order terms of Roemer delay. PEXO is able to avoid the bias introduced by decoupling the target system and the solar system and to account for the atmospheric effects that set a practical limit for ground-based radial velocities close to 1 cm/s. Considering the various caveats in barycentric correction and ancillary data required to realize cm/s modeling, we recommend the preservation of original observational data.