A new, highly optimized implementation of the Hiller-Sucher-Feinberg (HSF) identity is presented. The HSF identity, when applied to molecular wave functions calculated with Gaussian-type basis functions, not only improves the overall accuracy of the electron density by more than an order of magnitude, but also yields approximate cusps at nuclei. The three classes of molecular integrals, L, U, and V, which are encountered in the calculation of the HSF density, are derived in compact form. Efficient algorithms for the accurate evaluation of these integrals are detailed, including a novel approach to the necessary numerical quadratures and the thresholding of two-electron V integrals. Hartree-Fock (HF) electron densities calculated with both the conventional definition and from the HSF identity are compared to their respective HF limits for a variety of diatomic molecules and basis sets. The average error in the calculated HSF electron densities at non-hydrogen nuclei equals 0.17%, which constitutes a marked improvement over an error of 5.77% in the conventional densities.