Galerkin discretizations of integral equations in ${ℝ}^d$ require the evaluation of integrals $I={\int }_{S^{\left(1\right)}}{\int }_{S^{\left(2\right)}}g(x,y)\mathrm{d}y\mathrm{d}x$ where S(1),S(2) are d-simplices and g has a singularity at x = y. We assume that g is Gevrey smooth for x$\ne$y and satisfies bounds for the derivatives which allow algebraic singularities at x = y. This holds for kernel functions commonly occurring in integral equations. We construct a family of quadrature rules ${𝒬}_N$ using N function evaluations of g which achieves exponential convergence |I – ${𝒬}_N$| ≤C exp(–rNγ) with...