  We present a fast algorithm to calculate Coulomb/exchange integrals of
prolate spheroidal electronic orbitals, which are the exact solutions of the
single-electron, two-center Schr\"odinger equation for diatomic molecules. Our
approach employs Neumann's expansion of the Coulomb repulsion 1/|x-y|, solves
the resulting integrals symbolically in closed form and subsequently performs a
numeric Taylor expansion for efficiency. Thanks to the general form of the
integrals, the obtained coefficients are independent of the particular
wavefunctions and can thus be reused later.
  Key features of our algorithm include complete avoidance of numeric
integration, drafting of the individual steps as fast matrix operations and
high accuracy due to the exponential convergence of the expansions.
  Application to the diatomic molecules O2 and CO exemplifies the developed
methods, which can be relevant for a quantitative understanding of chemical
bonds in general.
