Singular integrals frequently arise in the calculation of the diagonal elements of BEM influence matrices and related multiple reciprocity applications. The underlying theory of Gaussian integration is applied to develop a general algorithm for the determination of the Gaussian parameters with the singular fundamental solution chosen as weight function. In contrast to singularities like ln(1/x), 1/x, and 1/x(2), which are frequently met in applied mechanics, the singularity of K-0(cx), the fundamental solution of the modified Helmholtz equation, is not treated extensively in the literature. The general algorithm is applied to the special case of this fundamental solution, and the Gaussian weights and ordinates are determined. Numerical experiments and comparisons with alternate methods of integration are carried out to assess the merit of the newly developed quadrature. (C) 2000 Elsevier Science Ltd. All rights reserved.