  A formalism for the numerical integration of one- and two-loop integrals is
presented. It is based on subtraction terms which remove the soft, collinear
and some of the ultraviolet divergences from the integrand. The numerical
integral is performed in the Feynman parameter space, using a complex contour
deformation to ensure robust convergence even in the presence of physical
thresholds. The application of the proposed procedure is demonstrated with
several one- and two-loop examples. An implementation in the program NICODEMOS
is publicly available, which currently incorporates only one-loop
functionality, but an extension to two-loop cases is planned for future
versions.
