A fully analytical implementation of the nuclear magnetic resonance (NMR) indirect nuclear spin-spin coupling constants at the density-functional theory (DFT) level is presented. The implementation involves all four contributions of the nonrelativistic Ramsey theory: The dia- and para-magnetic spin-orbit contributions as well as the paramagnetic Fermi-contact and spin-dipole contributions. Three different exchange-correlation functionals-LDA (local density approximation), BLYP (Becke-Lee-Yang-Parr), and B3LYP (hybrid BLYP)-are tested by comparison with experiment and high-level ab initio calculations for a series of molecules containing first-row elements. All three levels of theory represent a significant improvement on restrictred Hartree-Fock (RHF) theory in the sense that the RHF instability problems are absent in DFT. Also, there is a steady improvement in the quality of the calculated spin-spin couplings in the sequence LDA, BLYP, and B3LYP. For the first-row molecules investigated by us, the hybrid B3LYP functional performs particularly well, with errors similar to those observed at the best ab initio levels of theory.