  We present a new method for the numerical evaluation of arbitrary loop
integrals in dimensional regularization. We first derive Mellin-Barnes integral
representations and apply an algorithmic technique, based on the Cauchy
theorem, to extract the divergent parts in the epsilon->0 limit. We then
perform an epsilon-expansion and evaluate the integral coefficients of the
expansion numerically. The method yields stable results in physical kinematic
regions avoiding intricate analytic continuations. It can also be applied to
evaluate both scalar and tensor integrals without employing reduction methods.
We demonstrate our method with specific examples of infrared divergent
integrals with many kinematic scales, such as two-loop and three-loop box
integrals and tensor integrals of rank six for the one-loop hexagon topology.
