  The Dirac-Coulomb equation with positive-energy projection is solved using
explicitly correlated Gaussian functions. The algorithm and computational
procedure aims for a parts-per-billion convergence of the energy to provide a
starting point for further comparison and further developments in relation with
high-resolution atomic and molecular spectroscopy. Besides a detailed
discussion of the implementation of the fundamental spinor structure,
permutation and point-group symmetries, various options for the positive-energy
projection procedure are presented. The no-pair Dirac-Coulomb energy converged
to a parts-per-billion precision is compared with perturbative results for
atomic and molecular systems with small nuclear charge numbers. The subsequent
paper [Paper II: D. Ferenc, P. Jeszenszki, and E. M\'atyus (2022)] describes
the implementation of the Breit interaction in this framework.
