I'm curious to know how QUADPACK's QAG routine works. My understanding is that it begins by calculating on each subinterval the numerical quadrature with a Gaussian-Legendre rule and a nested Kronrod rule to provide estimates of the integral and the local error. Then, globally the interval with the highest error is bisected and the process is repeated. However, since Gauss-Legendre abscissae do not nest, can any of the previous function evaluations be reused at the next iteration? Or are all previous function evaluations discarded when an interval is bisected and reevaluated? If the latter is true this seems incredibly inefficient, so why would anyone choose such an adaptive quadrature scheme like that in QUADPACK over one where the function evaluations may be reused, such as Newton-Cotes or Clenshaw-Curtis?