The present work introduces a stable, semi-implicit, one stage integration scheme for rigid multibody systems subject to mixed holonomic and nonholonomic constraints as well as dry frictional impacts and contacts. A stable direct-iterative splitting scheme to solve the latter is also presented, and is shown to be suitable for real-time simulation of large multibody systems, such as those used for 3D graphics, simulation-based training systems. We use a Lagrangian framework in conjunction with the discrete-time D'Alembert's principle to introduce physically motivated singular perturbations which regularize and stabilize the numerical method. Lower bounds on the perturbations which guarantee numerical stability are provided, and their physical validity is demonstrated at the numerical level. The theoretical formulation uses massless ghost particles in the Lagrangians of mechanical systems. The coordinates and velocities of these converge strongly to Lagrange multipliers for holonomic and nonholonomic constraints, respectively, in the singular limit of zero relaxation. The ghost formulation allows for the systematic treatment of non-ideal, dissipative, constitutive laws such as Coulomb friction. A variational model for the latter is constructed and proven to be solvable in discrete time. Several splitting schemes are investigated mathematically and compared at the numerical level. Convergence and non-convergence properties of these are demonstrated mathematically and numerically.