This article summarizes the development of a fast boundary element method for the linear Poisson-Boltzmann equation governing biomolecular electrostatics. Unlike previous fast boundary element implementations, the present treatment accommodates finite salt concentrations thus enabling the study of biomolecular electrostatics under realistic physiological conditions. This is achieved by using multipole expansions specifically designed for the exponentially decaying Green's function of the linear Poisson-Boltzmann equation. The particular formulation adopted in the boundary element treatment directly affects the numerical conditioning and thus convergence behavior of the method. Therefore, the formulation and reasons for its choice are first presented. Next, the multipole approximation and its use in the context of a fast boundary element method are described together with the iteration method employed to extract the surface distributions. The method is then subjected to a series of computational tests involving a sphere with interior charges. The purpose of these tests is to assess accuracy and verify the anticipated computational performance trends. Finally, the salt dependence of electrostatic properties of several biomolecular systems (alanine dipeptide, barnase, barstar, and coiled coil tetramer) is examined with the method and the results are compared with finite difference Poisson-Boltzmann codes.