  We report a high-precision numerical estimation of the critical exponent
$\alpha$ of the specific heat of the random-field Ising model in four
dimensions. Our result $\alpha = 0.12(1)$ indicates a diverging specific-heat
behavior and is consistent with the estimation coming from the modified
hyperscaling relation using our estimate of $\theta$ via the anomalous
dimensions $\eta$ and $\bar{\eta}$. Our analysis benefited form a
high-statistics zero-temperature numerical simulation of the model for two
distributions of the random fields, namely a Gaussian and Poissonian
distribution, as well as recent advances in finite-size scaling and reweighting
methods for disordered systems. An original estimate of the critical slowing
down exponent $z$ of the maximum-flow algorithm used is also provided.
