The present contribution reviews the research over the past decade on the computation of the three-center nuclear attraction integrals over B functions, during which we have developed various numerical methods aimed at creating numerical algorithms capable of efficiently and accurately handling these pathological integrals. Given the enormous number of integrals involved, reducing calculation time while maintaining high precision is of paramount importance. One of the key advancements was the introduction of nonlinear transformations, particularly the D and G transformations, which have significantly outperformed previous methods, such as Gauss-Laguerre quadrature, Wynn’s epsilon algorithm and Levin’s u-transform, which are among the most popular convergence accelerators. These transformations D and G marked a substantial leap forward in the efficiency and accuracy of the calculations. Following this, by applying the S transformation, the spherical Bessel integral can be converted into a much more favorable sine integral or into a combination of sine and cosine integrals. This has led to methods that are not highly accurate but also the most efficient compared to existing alternatives. More recently, we introduced an approach that applies two types of double-exponential transformations to the resulting sine integral, producing a quadrature formula that is both highly accurate and efficient. This method enables the approximation of molecular integrals to a predefined high accuracy while significantly reducing computation time.