I think in the early days one really used just the normal integral schemes Obara-Saika, McMurchie-Davidson, Rys where for one exponent just a Gaussian with exponent zero was used. Later the people examined the schemes for this special purpose and modified them for calculating 3 index coulomb integrals. A quite recent paper on this topic can be found in Gyula Samu and Mihály Kállay, J. Chem. Phys. 2017, 146, DOI: 10.1063/1.4983393.
One minor thing I want to point out is that you presented RI/DF in a slightly unusual way (notation). Using RI/DF the integrals are evaluated as
$$
(ij|kl) \approx \sum_{PQ}^{N_{aux}} (ij|P)[V^{-1}]_{PQ} (Q|kl)
$$
where
$$
V_{PQ} = \left(P|Q\right) = \int \int {{\phi_{P}({r_1}) \frac{1}{r_{12}} \phi_{Q}({r_2})}}d{r_1}d{r_2}
$$
Of course you can rewrite it a bit with forming $V_{PQ}^{-1/2}$ to arrive at a similar expression as given by you.