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.