Bohr-Sommerfeld quantization provides an approximate recipe for recovering the spectrum of a quantum integrable system. Is there a mathematically rigorous explanation why this recipe works? In particular, I suppose it gives an exact description of the large quantum number asymptotics, which should be a theorem.
Also, is there a way to make the recipe more precise by adding corrections of some sort?
Yes, it can be made precise and corresponds to the leading order of the semiclassical expansion (WKB approximation) in $\hbar$. See Faddeev-Yakubovsky's "Lectures on quantum mechanics for mathematics students" (§20, formula (13)). An approach inspired by geometric quantization is explained in chapter 4 in Bates-Weinstein's Lectures on the Geometry of Quantization.
This answer addresses the geometrical origin of the Bohr-Sommerfeld condition. In geometric quantization, the additional structure required beyond the symplectic data of the phase space is a polarization. The quantization spaces are constructed as spaces of polarized sections with respect a polarization. The most "obvious" type of a polarization is the Kahler polarization, where the quantization spaces are spaces of holomorphic sections of a pre-quantum line bundle. Simple examples of systems which can be quantized by means of a Kahler polarization are the Harmonic oscillator and the spin. Another type of polarization is the real polarization (Please see for example the Blau's lecttures), which is locally equivalent to the polarization of a cotangent bundle. A real polarizations foilates the phase space (symplectic manifold) into Lagrangian submanifolds. When the leaves are compact, the quantum Hilbert space consists of sections with support only on certain leaves, which are exactly those which satisfy the Bohr-Sommerfeld condition. In this case, the quantum phase space is generated by distributional sections supported solely on the Bohr-Sommerfeld leaves (This result is due to Snyatycki). For example in the case of the spin, the Bohr-Sommerfeld leaves are small circles at half integer values of the $z$ coordinate in the two dimensional sphere. A more sophisticated example of Bohr-Sommerfeld leaves is the Gelfand-Cetlin system on flag manifolds.
Many classical phase admit both Kahler and real polarizations. It is interesting that in many cases the quantization Hilbert spaces are unitarily equivalent (i.e. the quantization is polarization independent). Please see for example Nohara's exposition.
Contrarily to what is generally believed, a semiclassical approximation is achieved through two different series: One is WKB series and the other is the Wigner-Kirkwood series, the latter being a gradient expansion. In both cases, eigenvalues are obtained by the Bohr-Sommerfeld rule but just at the leading order. I have proved this here (this paper appeared in Proceedings of Royal Society A). This proof is rigorous and quite different from what one finds on standard textbooks. Besides, it produces the full series for the exact eigenvalues with at leading order the ordinary Bohr-Sommerfeld rule.
Here's a very basic way to see this easily:
This is Action, search for the "Abbreviated Action", and it has the SI unit of Joule-second. This is the equation that Planck (and later, Einstein) used: $$E=nhf$$ for $n=1,2,3...$ and $f$ in frequency, in unit of 1/second).
This means that the Planck's constant has also an unit of Joule-second, therefore, you can interpret $nh$ as the Action of the system (since $n$ is dimensionless, it will maintain the unit, and it is only used to give the "correct" answer, since not every mechanical system has $h$ as the Action, $n$ should be used).
So $$nh=\int p\,dq$$ which is the Sommerfield rule for quantization. This is just an intuition of how it works.
perhaps it can be derived from te approximation over the density of states
$$ N(E)= \sum_{n=0}^{\infty}\theta(E-E_{n})\approx \frac{1}{2\pi \hbar}\iint_{V}\theta(E-H)dxdp $$
with $ H= P^{2}/2m +V(x) $ is the Hamiltonian of the particle and $ \theta (x) $ is heaviside step function
