We present efficient techniques to approximate singular and nearly singular surface integrals arising massively in Galerkin boundary element discretizations of Fredholm integral equations on two-dimensional surfaces. The technique is based on a new representation of the functionals which arise by computing the so-called local element matrices, i.e., integrals over pairs of panels. The remaining integrals are approximated by introducing relative co-ordinates which fix the location of the singularity. These integrals can be integrated using polar co-ordinates. It turns out that the number of kernel evaluations which is needed to compute the integrals up to the required accuracy is independent of the order of the singularity. This enables us to use the hypersingular formulation of integral equations which is the method of choice from the theoretical point of view, i.e., stability, robustness with respect to non-smooth surfaces, etc.