The computational methods for the diffraction integrals that occur in the Extended Nijboer-Zernike (ENZ-) approach to circular, aberrated, defocused optical systems are reviewed and updated. In the ENZ-approach, the Debye approximation of Rayleigh’s integral for the through-focus, complex, point-spread function is evaluated in semi-analytic form. To this end, the generalized pupil function, comprising phase aberrations as well as amplitude non-uniformities, is assumed to be expanded into a series of Zernike circle polynomials, and the contribution of each of these Zernike terms to the diffraction integral is expressed in the form of a rapidly converging series (containing power functions and/or Bessel functions of various kinds). The procedure of expressing the through-focus point-spread function in terms of Zernike expansion coefficients of the pupil function can be reversed and has led to the ENZ-method of retrieval of pupil functions from measured through-focus (intensity) point-spread functions. The review and update concern the computation for systems ranging from as basic as having low NA and small defocus parameter to high-NA systems, with vector fields and polarization, meant for imaging of extended objects into a multi-layered focal region. In the period 2002-2010, the evolution of the form of the diffraction integral (DI) was dictated by the agenda of the ENZ-team in which a next instance of the DI was handled by amending the computation scheme of the previous one. This has resulted into a variety of ad hoc measures, lack of transparency of the schemes, and sometimes prohibitively slow computer codes. It is the aim of the present paper to reconstruct the whole building of computation methods, using consistently more advanced mathematical tools. These tools are explicit Zernike expansion of the focal factor in the DI, Clebsch-Gordan coefficients for the omnipresent problem of linearizing products of Zernike circle polynomials, recursions for Bessel functions, binomials and for the coefficients of algebraic functions that occur as pre-factors of the focal factor in the DI. This results in a series representation of the DI involving (spherical) Bessel functions and Clebsch-Gordan coefficients, in which the dependence of the DI on parameters of the optical configuration, on focal values, on spatial variables in the image planes, and on degree and azimuthal order of the circle polynomials are separated. This separation of dependencies, together with bounds on Clebsch-Gordan coefficients and spherical Bessel functions, facilitate the error analysis for the truncation of series, showing that in the new scheme the DI can be computed virtually without loss-of-digits. Furthermore, this separation allows for a modular implementation of the computation scheme that offers speed and flexibility when varying the various parameters and variables. The resulting scheme is pre-eminently appropriate for use in advanced optical simulations, where large defocus values, high NA and Zernike terms of high order and degree occur.