In rigid multibody dynamics simulation using absolute coordinates, a choice must be made in relation to how to keep track of the attitude of a body in 3D motion. The commonly used choices of Euler angles and Euler parameters each have drawbacks, e.g., singularities, and carrying along extra normalization constraint equations, respectively. This contribution revisits an approach that works directly with the orientation matrix and thus eschews the need for generalized coordinates used at each time step to produce the orientation matrix A. The approach is informed by the fact that rotation matrices belong to the SO(3) Lie matrix group. The numerical solution of the dynamics problem is anchored by an implicit first order integration method that discretizes, without index reduction, the index 3 Differential Algebraic Equations (DAEs) of multibody dynamics. The approach handles closed loops and arbitrary collections of joints. Our main contribution is the outlining of a systematic way for computing the first order variations of both the constraint equations and the reaction forces associated with arbitrary joints. These first order variations in turn anchor a Newton method that is used to solve both the Kinematics and Dynamics problems. The salient observation is that one can express the first order variation of kinematic quantities that enter the kinematic constraint equations, constraint forces, external forces, etc., in terms of Euler infinitesimal rotation vectors. This opens the door to a systematic approach to formulating a Newton method that provides at each iteration an orthonormal rotation matrix A. The Newton step calls for repeatedly solving linear systems of the form Gδ = e, yet evaluating the iteration matrix G and residuals e is inexpensive, to the point where in the Part 2 companion contribution the proposed formulation is shown to be two times faster for Kinematics and Dynamics analysis when compared to the Euler parameter and Euler angle approaches in conjunction with a set of four mechanisms.