We describe an efficient biased Monte-Carlo method for calculating the diag
rams appearing in the coefficients of the so-called bridge function B = Sig
ma(n=2)(infinity) b(n)rho(n) of the integral equation theory of liquids. Th
ese diagrams represent multi-dimensional integrals of products of 'bond' fu
nctions of the intermolecular distances. The method rests on the generation
of independent Markov chains and is well adapted to highly parallel comput
ation. It can be used for systems with any pair potential. The feasibility
and efficiency of the method are demonstrated for the second and third orde
r coefficients of the bridge functions of fluids of hard and Lennard-Jones
spheres. For these systems there are analytical expressions of the bridge f
unction deduced from computer simulations to which we compare our bridge fu
nction approximations which include the second and third order coefficients
with h as the bond function. Our new approximations of the bridge function
are used in the closure of the Ornstein-Zernike relation. The obtained str
uctural and thermodynamical properties are found in better agreement with t
he exact simulation data than the hypernetted chain results.