  We present a precise and efficient computation of the two-loop amplitudes
entering the Higgs boson pair production process via gluon fusion. Our approach
is based on the small-Higgs-mass expansion while keeping the full dependence on
the top quark mass and other kinematic invariants. We compare our results to
the up-to-date predictions based on a combination of sector decomposition and
high-energy expansion. We find that our method provides precision numeric
predictions in the entire phase space, while at the same time is highly
efficient as the computation can be easily performed on a normal desktop or
laptop computer. Our method is valuable for practical phenomenological studies
of the Higgs boson pair production process, and can also be applied to other
similar processes.
