Boron nitride (BN) is a material with outstanding technological promise because of its exceptional thermochemical stability, structural, electronic and thermal conductivity properties, and extreme hardness. Yet, the relative thermodynamic stability of its most common polymorphs (diamond-like cubic and graphite-like hexagonal) has not been resolved satisfactorily because of the crucial role played by kinetic factors in the formation of BN phases at high temperatures and pressures (experiments), and by competing bonding, electrostatic and many-body dispersion forces in BN cohesion (theory). This lack of understanding hampers the development of potential technological applications, and challenges the boundaries of fundamental science. Here, we use high-level first-principles theories that correctly reproduce all important electronic interactions (the adiabatic-connection fluctuation-dissipation theorem in the random phase approximation) to estimate with unprecedented accuracy the energy differences between BN polymorphs, and thus overcome the accuracy hurdle that hindered previous theoretical studies. We show that the ground-state phase of BN is cubic and that the frequently observed two-dimensional hexagonal polymorph becomes entropically stabilized over the cubic at temperatures slightly above ambient conditions (Tc = 63+-20'C). We also reveal a new low-symmetry monoclinic phase that is extremely competitive with the other low-energy polymorphs and which could explain the origins of the experimentally observed ``compressed h--BN'' phase. Our theoretical findings therefore should stimulate new experimental efforts in bulk BN as well as promote the use of high-level theories in modelling of technologically relevant van der Waals materials.<br>