We provide a fully numerical, deterministic integration at the level of the three- and four-point functions, in the reduction of the one-loop hexagon integral by sector decomposition. For the corresponding two- and three-dimensional integrals we use an adaptive numerical approach applied recursively in two and three dimensions, respectively. The adaptive integration is coupled with an extrapolation for an accurate, automatic treatment of integrand singularities arising from vanishing denominators in the interior of the integration domain. Furthermore, the recursive procedure alleviates extensive memory use as incurred with standard adaptive, multidimensional integration software. Tensor integrals are handled automatically by this technique and the separation of infrared singularities follows naturally by dimensional regularization.