How to calculate the overlap of four cartesian gaussian function i.e $\int \phi_{A}(r)\phi_{B}(r)\phi_{C}(r)\phi_{D}(r)d^{3}r$?

Question

I want to calculate $\int \phi_{A}(r)\phi_{B}(r)\phi_{C}(r)\phi_{D}(r)d^{3}r$ where $\phi_{A}(r)=(x-A_{x})^{l1}(y-A_{y})^{l2}(z-A_{z})^{l3}e^{-\alpha |r-A|^{2}}$. l1, l2 and l3 are integers Can you suggest any book and any code in c++/matlab that can calculate this quantity?