We present a method for solving impurity models with electron-phonon coupling, which treats the phonons efficiently and without approximations. The algorithm is applied to the Holstein-Hubbard model in the dynamical mean field approximation, where it allows access to strong interactions, very low temperatures, and arbitrary fillings. We show that a renormalized Migdal-Eliashberg theory provides a reasonable description of the phonon contribution to the electronic self-energy in strongly doped systems, but fails if the quasiparticle energy becomes of order of the phonon frequency.