  One-loop two-, three- and four-point scalar functions are analytically
integrated directly such that they are expressed in terms of Lauricella's
hypergeometric function $F_D$. For two- and three-point functions, exact
expressions are obtained with arbitrary combination of kinematic and mass
parameters in arbitrary space-time dimension. Four-point function is expressed
in terms of $F_D$ up to the finite part in the expansion around 4-dimensional
space-time with arbitrary combination of kinematic and mass parameters. Since
the location of the possible singularities of $F_D$ is known, information about
the stabilities in the numerical calculation is obtained. We have developed a
numerical library calculating $F_D$ around 4-dimensional space-time. The
numerical values for IR divergent cases of four-point functions in massless QCD
are calculated and agreed with golem95 package.

High Energy Physics - Phenomenology