I would like to numerically integrate a smooth function (to be precise, I am computing differential entropies of Gaussian mixtures of the form $\int f(x) \log(1/f(x)) dx$).

Is there any software package or library that gives me the exact answer up to X digits (and potentially more), where X can be rather large? The GSL numerical integration seems to provide a solution with its claim:

The algorithms attempt to estimate the absolute error ABSERR = |RESULT - I| in such a way that the following inequality holds,

|RESULT - I| <= ABSERR <= max(epsabs, epsrel |I|)

but first I don't understand what is exactly meant by "attempt" (boldface added above) and second I would like to go beyond machine floating point if possible.

I am OK if the software gives up (if the requested accuracy/precision is too severe), but once an answer is out, it ought to be rather reliable. (Otherwise, what would be the point of computing it to begin with?)

This must be a very basic problem in numerical integration, but my search over the past few weeks have been futile. Any thread would be greatly appreciated.

You can't bound the accuracy of any numerical integration scheme without some control over how crazy the integrand can be. Smoothness is obviously helpful, but if the derivatives aren't bounded, then it's possible to introduce a $C^{\infty}$ bump into the integrand that can have an arbitrarily large impact on the value of the integral and you might not ever sample from that portion of the integration interval.

In practice, adaptive quadrature schemes attempt to estimate the error and these error estimates are quite reliable for most practical situations in which the integrand is smooth and the derivatives are bounded (even if you don't know the bounds), but there are no theoretical guarantees for a general smooth integrand.

I would suggest that you look at symbolic computation software packages like Maple and Mathematica that offer arbitrarily high precision floating point arithmetic and have robust numerical quadrature routines.