I am interested in such extension.
Do you plan to make it precision-independent, so that it can compute with any desired precision using multiple-precision scalar types (e.g. "mpreal")?
(most of the code in Eigen is precision-independent and "mpreal" is included in distribution).
The main obstacle for this is that QUADPACK uses pre-computed coefficients in "double" precision only.
So that we have to implement additional algorithms to compute Gauss-Kronrod nodes & weights for (any) required precision on the fly.
The are several algorithms to generate GK:
D. Laurie (1997). Calculation of Gauss-Kronrod Quadrature Rules. Mathematics of Computation, 66(219).
G. Monegato (1978).. Some remarks on the construction of extended Gaussian quadrature rules. Mathematics of Computation, 32(141).
Here is one of the (unfinished) attempts to implement Monegato method:
I am ready to help with any "mpreal"-related questions (as I am its author).