  In this article we propose a new adaptive numerical quadrature procedure
which includes both local subdivision of the integration domain, as well as
local variation of the number of quadrature points employed on each
subinterval. In this way we aim to account for local smoothness properties of
the function to be integrated as effectively as possible, and thereby achieve
highly accurate results in a very efficient manner. Indeed, this idea
originates from so-called hp-version finite element methods which are known to
deliver high-order convergence rates, even for nonsmooth functions.
