  A scheme for systematically achieving accurate numerical evaluation of
multi-loop Feynman diagrams is developed. This shows the feasibility of a
project aimed to produce a complete calculation for two-loop predictions in the
Standard Model. As a first step an algorithm, proposed by F. V. Tkachov and
based on the so-called generalized Bernstein functional relation, is applied to
one-loop multi-leg diagrams with particular emphasis to the presence of
infrared singularities, to the problem of tensorial reduction and to the
classification of all singularities of a given diagram. Successively, the
extension of the algorithm to two-loop diagrams is examined. The proposed
solution consists in applying the functional relation to the one-loop
sub-diagram which has the largest number of internal lines. In this way the
integrand can be made smooth, a part from a factor which is a polynomial in the
vector of Feynman parameters needed for the complementary sub-diagram with the
smallest number of internal lines. Since the procedure does not introduce new
singularities one can distort the integration hyper-contour into the complex
hyper-plane, thus achieving numerical stability. The algorithm is then modified
to deal with numerical evaluation around normal thresholds. Concise and
practical formulas are assembled and presented, numerical results and
comparisons with the available literature are shown and discussed for the
so-called sunset topology.
