A novel approach to the description of superconductors in thermal equilibrium is developed within a formally exact density-functional framework. The theory is formulated in terms of three ``densities'': the ordinary electron density, the superconducting order parameter, and the diagonal of the nuclear N-body density matrix. The electron density and the order parameter are determined by Kohn-Sham equations that resemble the Bogoliubov-de Gennes equations. The nuclear density matrix follows from a Schroedinger equation with an effective N-body interaction. These equations are coupled to each other via exchange-correlation potentials which are universal functionals of the three densities. Approximations of these exchange-correlation functionals are derived using the diagrammatic techniques of many-body perturbation theory. The bare Coulomb repulsion between the electrons and the electron-phonon interaction enter this perturbative treatment on the same footing. In this way, a truly ab-initio description is achieved which does not contain any empirical parameters.