This paper presents an equation of motion for numerical simulation of constrained mechanical systems with holonomic and nonholonomic constraints. In order to avoid the error accumulation typically experienced in such simulations, the standard equation of motion is enhanced with embedded force and impulse terms which perform continuous constraint and energy correction along the numerical solution. To avoid interference between the kinematic constraint correction and the energy correction terms, both are derived by taking the geometry of the constrained dynamics rigorously into account. In this light, enforcement of the (ideal) holonomic and nonholonomic kinematic constraints are performed using ideal forces and impulses, while the energy conservation law is considered as a nonideal nonlinear nonholonomic constraint on the simulated motion, and as such it is enforced with nonideal forces. As derived, the equation can be directly discretized and integrated with an explicit ODE solver avoiding the need for expensive implicit integration and iterative constraint stabilization. Application of the proposed equation is demonstrated on a representative example. A more elaborate discussion of practical implementation is presented in Part II of this work.