Two-electron integral algorithm

I am writing a molecular integral code from scratch. Right now my code have the following structure:

loop over shellpair ab:
    loop over shellpair cd:
        loop over primitive ij in ab:
            calculate E_ij
            loop over primitive kl in cd:
                calculate E_kl
                calculate R_ij_kl
                calculate V_ab_cd_ij_kl
        contract V_ab_cd_ij_kl to V_ab_cd

The above pseudo-code have alot of implicit notation, so let me explain what it means before asking my question.

E_ij and E_kl is the expansions coefficients from the McMurchie Davidson scheme.
R_ij_kl is the Hermite integrals also used in the Mcmurchie Davidson scheme.
V_ab_cd_ij_kl is the cartesian integrals found by (note the equations below are not strictly correct):

$$V_{ab,cd,ij,kl}=sum E_{kl}sum E_{ij} R_{ij,kl}$$

And V_ab_cd is the contraction of the primitive to the basisfunction:

$$V_{ab,cd}=sum_isum_jsum_ksum_l c_ic_jc_kc_lV_{ab,cd,ij,kl}$$

When calculating integrals of higher angular momentum (above $0$), I make sure to reuse common elements in E and R. I.e. when doing p-type integrals px, py and pz have the first element in common, and it is only calculated once in the recurrence relations of the McMurchie Davidson scheme.

For my implementation of the above pseudo-code even the (ss|ss) integrals are very slow. When doing benenze with 6-31G the (ss|ss) integrals uses 2.5 seconds, which is more than can be justified by the use of Python (accelerated with Numba) instead of a compiled language. I have thought about rearranging and saving more intermidiates in memory, i.e. doing something like:

 loop over shellpair ab:
    calculate E_ab (Here E_ab means all E_ij saved in memory)
    loop over shellpair cd:
        calculate E_cd (Here E_ab means all E_kl saved in memory)
        loop over primitive ij in ab:
            loop over primitive kl in cd:
                calculate R_ij_kl
                calculate V_ab_cd_ij_kl
        contract V_ab_cd_ij_kl to V_ab_cd

In the above by calculating E_ab and E_cd, recomputation of E_ij and E_kl is avoided. At this point I know that I lag some expertise in knowing what is worth it compared to memory usage.

My question is know, what is a proper way of implementing molecular integrals with regard of what intermediates to keep at which points. (I know thing like vertical recurrence relations exist, but these will for obvious reasons not speed on my (ss|ss) integrals).