  Seniority-zero geminal wavefunctions are known to capture bond-breaking
correlation. Among this class of wavefunctions, Richardson-Gaudin states stand
out as they are eigenvectors of a model Hamiltonian. This provides a clear
physical picture, clean expressions for reduced density matrix (RDM) elements,
and systematic improvement (with a complete set of eigenvectors). Known
expressions for the RDM elements require the computation of rapidities, which
are obtained by first solving for the so-called eigenvalue based variables
(EBV) then root-finding of a Lagrange interpolation polynomial. In this
manuscript we obtain expressions for the RDM elements directly in terms of the
EBV. The final expressions can be computed with the same cost as the rapidity
expressions. Therefore, except in particular circumstances, it is entirely
unnecessary to compute rapidities at all. The RDM elements require numerically
inverting a matrix and while this is usually undesirable we demonstrate that it
is stable, except when there is degeneracy in the single-particle energies. In
such cases a different construction would be required.
