As it turns out, there is a better approach using Clenshaw's algorithm! A review of different methods is given in arXiv:1410.1748 [physics.chem-ph]. Now I could get up to high values of \(\ell\) and \(m\) without everything blowing up! From the ALPs we can then simply obtain the SHs as
\[
Y_\ell^m(\theta, \phi) = \sqrt{\frac{(2 \ell + 1) (\ell - m)!}{4\pi (\ell + m)!}} P_\ell^m (\cos \theta) e^{i m \phi}.
\]
The final code can be viewed at https://www.shadertoy.com/view/3dKGWG, which borrows from iq's shader
@tpfto Thanks! Your posts have helped me many times over the years
@roninkoi Very nicely done! I had planned to write an improvement of Quilez's original demo after writing my own implementation of Clenshaw for other orthogonal polynomials (e.g. https://www.shadertoy.com/view/tlX3D7), but real life got in the way. Thanks for making it so I no longer have to. ;)