  We consider QCD radiative corrections to Higgs boson pair production through
gluon fusion in proton collisions. We combine the exact next-to-leading order
(NLO) contribution, which features two-loop virtual amplitudes with the full
dependence on the top quark mass $M_t$, with the next-to-next-to-leading order
(NNLO) corrections computed in the large-$M_t$ approximation. The latter are
improved with different reweighting techniques in order to account for
finite-$M_t$ effects beyond NLO. Our reference NNLO result is obtained by
combining one-loop double-real corrections with full $M_t$ dependence with
suitably reweighted real--virtual and double-virtual contributions evaluated in
the large-$M_t$ approximation. We present predictions for inclusive cross
sections in $pp$ collisions at $\sqrt{s}$=13, 14, 27 and 100TeV and we discuss
their uncertainties due to missing $M_t$ effects. Our approximated NNLO
corrections increase the NLO result by an amount ranging from +12% at
$\sqrt{s}$=13TeV to +7% at $\sqrt{s}$=100TeV, and the residual uncertainty from
missing $M_t$ effects is estimated to be at the few percent level. Our
calculation is fully differential in the Higgs boson pair and the associated
jet activity: we also present predictions for various differential
distributions at $\sqrt{s}$=14 and 100TeV. Our results represent the most
advanced perturbative prediction available to date for this process.
