We present the calculation of the cross section and invariant mass distribution for Higgs boson pair production in gluon fusion at next-to-leading order (NLO) in QCD. Top-quark masses are fully taken into account throughout the calculation. The virtual two-loop amplitude has been generated using an extension of the program GoSam supplemented with an interface to Reduze for the integral reduction. The occurring integrals have been calculated numerically using the program SecDec. Our results, including the full top-quark mass dependence for the first time, allow us to assess the validity of various approximations proposed in the literature, which we also recalculate. We find substantial deviations between the NLO result and the different approximations, which emphasizes the importance of including the full top-quark mass dependence at NLO.