How does one compute the number of unique 2-electron integrals for a given basis set?

I have run into trouble with an exercise in Szabo and Ostlund's Modern Quantum Chemistry.

Exercise 3.14: Assume the basis functions are real and use the symmetry of the two-electron integrals [$(\mu\nu|\lambda\sigma)=(\nu\mu|\lambda\sigma)=(\nu\mu|\sigma\lambda)$, etc.] to show that for a basis set of size $K=100$ there are $12,753,775=O(K^4/8)$ unique two-electron integrals.

I have seen answers to this on the Physics SE here and here, but the case-by-case analysis used in the answers seems confusing to me and not easily generalized to other values of $K$. The symmetry of the integrals when all the basis functions are real suggests that there should be a simpler way of determining this. How can one find the number of unique two-electron integrals in this case and, ideally, in general?

To find the number of unique $2e^-$ integrals $\left<AB|CD\right>$, it is useful to first find the number of unique $1e^-$ integrals $\left<A|B\right>$.

If $n$ is the number of real basis functions, we can find the number of $1e^-$ integrals as ${n+1}\choose{2}$. This stems from the general formula for combinations with replacement: $${n+k-1}\choose{k}$$ where $n$ is the selection pool and $k$ is how many are being chosen. Combinations with replacement are needed to ensure we count the cases where the two Gaussians are the same.

As an aside, this is only true for real basis functions, since we can only say the order doesn't matter (what we do when we compute the number of combinations) because all the basis functions are explicitly real and thus the complex conjugate equals the original function. This explains one of my own questions on the site regarding why we don't commonly use complex basis functions.

For this particular case ($n=100$), we find that the number of unique $1e^-$ integrals is $${{101}\choose{2}}=\frac{(101)*(100)}{2}=5050$$

So how does this allow us to find the number of unique $2e^-$ integrals? Well, for the $1e^-$ integrals, all we did was find the number of unique pairs of real basis functions. So, for a $2e^-$ integral, we might expect that all we need to do is find the unique pairs of basis function pairs. This winds up being exactly the case (again, only for real basis functions).

To find the pairs of pairs, we need to again find the number of combinations with replacement, but instead of sampling from the $100$ basis functions, we are now sampling from the $5050$ basis function pairs. We again need replacement to account for the fact that we can pick the same pair twice. So, we obtain: $${{5051}\choose {2}} =\frac{(5051)*(5050)}{2}=12,753,775$$ as the number of unique $2e^-$ integrals. This method avoids the checking of specific cases (i.e none of the basis functions are the same, two are the same, etc.) and gives a general way of finding the number unique integrals for any amount of real basis functions.

To show explicitly why complex basis functions aren't used, we can look at the same problem, but assume all the functions are complex instead. In bra-ket notation, the bra is the complex conjugate of the function inside, so for $\left<A|B\right>$, the order would matter if $A$ and $B$ are complex. So finding the number of pairs requires finding the number of permutations rather than combinations. This only requires a subtle change of the original calculation $$2!{101\choose2}= (101)*(100)=10100$$ So, we now have twice as many unique pairs, so there are twice as many $1e^-$ integrals. The equation for the $2e^-$ integrals remains the same, but with the new amount of unique pairs filled in. $${{10101}\choose {2}} =\frac{(10101)*(10100)}{2}=51,010,050$$ We still need to do a combination for the pairs because regardless of whether the functions are complex or real $\left<AB|CD\right>=\left<CD|AB\right>$ via exchange of electron labels. So using complex basis functions would double the number of $1e^-$ integrals and nearly quadruple the number of $2e^-$ integrals. So even if some feature of complex basis functions made them easier to integrate, it would have to make them easier enough to warrant doing and storing 2-4 times as many integrals.