  We present a calculation of the next-to-leading order QCD corrections to the
production of Higgs boson pairs in gluon fusion keeping the full dependence on
the mass of the top quark. The virtual corrections, involving two-loop
integrals with up to four mass scales, have been calculated numerically and we
present an efficient algorithm to obtain accurate results of the virtual
amplitude using numerical integrations. Taking the top quark mass into account
we obtain significant differences compared to results obtained in the heavy top
limit.
