We present a first-principles scheme for incorporating many-body interactions into the unified description of the quadratic optical response to light of noncentrosymmetric crystals. The proposed method is based on time-dependent current-density response theory and includes the electron-hole attraction via a tensorial long-range exchange-correlation kernel, which we calculate using the parameter-free bootstrap approximation. By bridging with the Wannier-interpolation of the independent-particle transition matrix elements, the resulting numerical scheme is very general and allows us to resolve narrow many-body spectral features at low computational cost. We showcase its potential by inspecting the second-harmonic generation in the benchmark zinc-blende semiconductor GaAs, the layered graphitic semiconductor $\mathrm{B}{\mathrm{C}}_{2}\mathrm{N}$, and the Weyl semimetal TaAs. Our results show that excitonic effects can give rise to large and sharply localized one- and two-photon resonances that are absent in the independent-particle approximation. We find overall good agreement with available experimental measurements, capturing the magnitude and peak structure of the spectrum as well as the angular dependence at fixed photon energy. The implementation of the method in Wannier-based code packages can serve as a basis for performing accurate theoretical predictions of quadratic optical properties in a vast pool of materials.