The thermodynamic properties of the Shastry-Sutherland model have posed one of the longest-lasting conundrums in frustrated quantum magnetism. Over a wide range on both sides of the quantum phase transition (QPT) from the dimer-product state to the plaquette-based ground state, neither analytical nor any available numerical methods have come close to reproducing the physics of the excited states and thermal response. We solve this problem in the dimer-product phase by introducing two qualitative advances in computational physics. One is the use of thermal pure quantum (TPQ) states to augment dramatically the size of clusters amenable to exact diagonalization. The second is the use of tensor-network methods, in the form of infinite projected entangled-pair states (iPEPS), for the calculation of finite-temperature quantities. We demonstrate convergence as a function of system size in TPQ calculations and of bond dimension in our iPEPS results, with complete mutual agreement even extremely close to the QPT. Our methods reveal a remarkably sharp and low-lying feature in the magnetic specific heat, whose origin appears to lie in a proliferation of excitations composed of two-triplon bound states. The surprisingly low energy scale and apparently extended spatial nature of these states explain the failure of less refined numerical approaches to capture their physics. Both of our methods will have broad and immediate application in addressing the thermodynamic response of a wide range of highly frustrated magnetic models and materials.