Many software packages provide means to calculate the Coulomb and exchange integrals, for example, the exchange
$$\int \phi^*(1) \, \psi^*(2) \frac{1}{|\mathbf{r_1} - \mathbf{r_2}|} \, \psi(1) \, \phi(2) \, d\mathbf{r_1} \, d\mathbf{r_2} $$
using gaussian type orbitals, either in cartesian or spherical coordinates.
I would like to know how this integral is evaluated in practice, since the Coulomb term $1/|\mathbf{r_1}-\mathbf{r_2}|$ should diverge whenever two different electron densities' volume elements overlap. I understand that when you use spherical coordinates, you can get rid of $1/r_{12}$ because the spherical volume element becomes something like $r^2 \cos{\theta} \, dr \, d\theta \, d\phi$, so you cancel $r^2$ in the numerator with $r$ in the denominator. But I'm having trouble to see how can this be done for cartesian GTO in cartesian coordinates.