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.