Background: The nuclear many-body system is a strongly correlated quantum system, posing serious challenges for perturbative approaches starting from uncorrelated reference states. The last decade has witnessed considerable progress in the accurate treatment of pairing correlations, one of the major components in medium-sized nuclei, reaching accuracies below the 1% level of the correlation energy. Purpose: Development of a quantum many-body method for pairing correlations that is (a) competitive in the 1% error range, and (b) can be systematically improved with a fast (exponential) convergence rate. Method: The present paper capitalizes upon ideas from Richardson-Gaudin integrability. The proposed method is a two-step approach. The first step consists of the optimization of a Richardson-Gaudin ground state as variational trial state. At the second step, the complete set of excited states on top of this Richardson-Gaudin ground state is used as an optimal basis for a Configuration Interaction method in an increasingly large effective Hilbert space. Results: The performance of the variational Richardson-Gaudin (varRG) and Richardson-Gaudin Configuration Interaction (RGCI) method is benchmarked against exact results using an effective $G$-matrix interaction for the Sn region. The varRG already reaches accuracies around the 1% level of the correlation energies, and the RGCI step sees an additional improvement scaling exponentially with the size of the effective Hilbert space. Conclusions: The Richardson-Gaudin models of integrability provide an optimized complete basis set for pairing correlations.