An alternative basis set for analytical computation of two electron integrals
A few things come to mind: 1) my understanding is that (ab|cd) in a gaussian basis is exactly solved for arbitrary "angular momenta" using recurrence relations such as the Obara-Saika method. The O-S method also has ways of generating the requisite derivatives with respect to coefficients, without which it is difficult to minimize the energy or do anything practical.

2) On the question of alternative basis sets to the Gaussian, one simply needs a complete basis set (of which there are many), helps to be orthogonal (almost as many), has to describe the boundary conditions of vanishing at infinity or being periodic with unit cell (fewer), hopefully is cusp-like near the nuclei (very few), and integration has to be tractable without resorting to quadrature, since there are four indices (each with 3 dimensions) but also millions of these integrals to compute -- and that leaves you with the Gaussian.