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 Mt, with the next-to-next-to-leading order (NNLO) corrections computed in the large-Mt approximation. The latter are improved with different reweighting techniques in order to account for finite-Mt effects beyond NLO. Our reference NNLO result is obtained by combining one-loop double-real corrections with full Mt dependence with suitably reweighted real--virtual and double-virtual contributions evaluated in the large-Mt approximation. We present predictions for inclusive cross sections in pp collisions at s√=13, 14, 27 and 100TeV and we discuss their uncertainties due to missing Mt effects. Our approximated NNLO corrections increase the NLO result by an amount ranging from +12% at s√=13TeV to +7% at s√=100TeV, and the residual uncertainty from missing Mt 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 s√=14 and 100TeV. Our results represent the most advanced perturbative prediction available to date for this process.