We address the fundamental problem of computing range functions for a real function f: R→ R. In our previous work , we introduced recursive interpolation range functions based on the Cornelius–Lohner (CL) framework of decomposing f as f= g+ R, which requires to compute g(I) “exactly” for an interval I. There are two problems: this approach limits the order of convergence to 6 in practice, and exact computation is impossible to achieve in standard implementation models. We generalize the CL framework by allowing g(I) to be approximated by strong range functions, where ε> 0 is a user-specified bound on the error. This new framework allows, for the first time, the design of interval forms for f with any desired order of convergence. To achieve our strong range functions, we generalize Neumaier’s theory of constructing range functions from expressions over a Lipschitz class Ω of primitive functions. We show that the class Ω is very extensive and includes all common hypergeometric functions. Traditional complexity analysis of range functions is based on individual evaluation on an interval. Such analysis cannot differentiate between our novel recursive range functions and classical Taylor-type range functions. Empirically, our recursive functions are superior in the “holistic” context of the root isolation algorithm Eval. We now formalize this holistic approach by defining the amortized complexity of range functions over a subdivision tree. Our theoretical model agrees remarkably well with the empirical results. Among our previous novel range functions, we identified a Lagrange-type range function as the overall winner. In this paper, we introduce a Hermite-type range function that is even better. We further explore speeding up applications by choosing non-maximal recursion levels.