The boundary element method requires accurate calculation of line and surface integrals in which a factor such as ln (1/r), 1/r or 1/r2 is ill-behaved when r is nearly zero over a small portion of the domain. This paper shows that the underlying theory of Gaussian integration can be applied to get a quadrature formula for any particular distance from the element. An accurate integral can then be calculated with very few evaluation points by (a) computing the distance to the element; (b) interpolating in a table to obtain quadrature points and weights; and (c) summing the weighted values of the non-singular part of the integrand. This provides a fast, accurate value of the integral with no uncertainty as to the number of evaluation points that are required. Line integrals with near-singular 1/r, 1/r2, and ln(1/r) terms are used as examples. Accuracy of the formulas is verified experimentally for a wide range of distances. Interpolation in tabulated data is shown to introduce only minor errors. The line integral formulas are applied to a triangular region, yielding exact formulas of any desired degree for the case where a 1/r singularity is at a vertex of the triangle.