We have carried out large-scale computational quantum chemistry calculations on the K computer to obtain heats of formation for C-60 and some higher fullerenes with the DSD-PBE-PBE/cc-pVQZ double-hybrid density functional theory method. Our best estimated values are 2520.0 +/- 20.7 (C-60), 2683.4 +/- 17.7 (C-70), 2862.0 +/- 18.5 (C-76), 2878.8 +/- 13.3 (C-78), 2946.4 +/- 14.5 (C-84), 3067.3 +/- 15.4 (C-90), 3156.6 +/- 16.2 (C-96), 3967.7 +/- 33.4 (C-180), 4364 (C-240) and 5415 (C-320) kJ mol(-1). In our assessment, we also find that the B3-PW91-D3BJ and BMK-D3(BJ) functionals perform reasonably well. Using the convergence behavior for the calculated per-atom heats of formation, we obtained the formula Delta H-f per carbon = 722n(-0.72) + 5.2 kJ mol(-1) (n = the number of carbon atoms), which enables an estimation of Delta H-f for higher fullerenes more generally. A slow convergence to the graphene limit is observed, which we attribute to the relatively small proportion of fullerene carbons that are in "low-strain" regions. We further propose that it would take tens, if not hundreds, of thousands of carbons for a fullerene to roughly approach the limit. Such a distinction, may be a contributing factor to the discrete properties between the two types of nanomaterials. During the course of our study, we also observe a fairly reliable means for the theoretical calculation of heats of formation for medium-sized fullerenes. This involves the use of isodesmic-type reactions with fullerenes of similar sizes to provide a good balance of the chemistry and to minimize the use of accompanying species.