Numerical integration is an operation that is frequently available in multiple precision numerical software packages. The different quadrature schemes used are considered well studied but the rounding errors that result from the computation are often neglected, and the actual accuracy of the results are therefore seldom rigorously proven.

We propose an implementation of the Gauss-Legendre quadrature scheme with bounded error: given a bound on the derivatives of a function we are able to compute an interval containing the true value of the integral, in arbitrary precision. The error analysis is given as well as experimental error measurements and timings, and a complete quadrature example.