Complexity in computing the integrals involving higher angular momentum I am quite familiar with Hartree Fock theory. However, writing it into program is an entirely different task, and that is where I have a couple of questions.

I wrote a Hartree-Fock routine in Python a couple of years ago. It involved a ultra-minimal basis: only one gaussian for two hydrogen atoms each. Kinetic energy term was pretty simple, so were the one-electron terms. Coulomb and exchange integrals were tricky, which involved evaluation of an erf. Now, my question is, if I included higher angular momentum basis functions (p,d,f,g type functions), what happens? Is it difficult to evaluate the J (Coulomb) and K (exchange) integrals because of the mathematical complexity? Or something else? Can someone give the math behind this and explain why it is difficult to evaluate?
It appears that the integrals that are difficult are the two-electron repulsion integrals. There exists a recursive relationship that can turn integrals with high-angular momentum functions into integrals with only s-type functions. Refer to this lecture for details.

See the following articles for details concerning these recurrence relations:

M. Head Gordon and J. Pople, J. Chem. Phys. 89, 5777 (1988)

S. Obara and A. Saika, J. Chem. Phys. 84, 3963 (1986)