  Computation of the spherical harmonic rotation coefficients or elements of
Wigner's d-matrix is important in a number of quantum mechanics and
mathematical physics applications. Particularly, this is important for the Fast
Multipole Methods in three dimensions for the Helmholtz, Laplace and related
equations, if rotation-based decomposition of translation operators are used.
In these and related problems related to representation of functions on a
sphere via spherical harmonic expansions computation of the rotation
coefficients of large degree $n$ (of the order of thousands and more) may be
necessary. Existing algorithms for their computation, based on recursions, are
usually unstable, and do not extend to $n$. We develop a new recursion and
study its behavior for large degrees, via computational and asymptotic
analyses. Stability of this recursion was studied based on a novel application
of the Courant-Friedrichs-Lewy condition and the von Neumann method for
stability of finite-difference schemes for solution of PDEs. A recursive
algorithm of minimal complexity $O\left(n^{2}\right)$ for degree $n$ and
FFT-based algorithms of complexity $O\left(n^{2}\log n\right) $ suitable for
computation of rotation coefficients of large degrees are proposed, studied
numerically, and cross-validated. It is shown that the latter algorithm can be
used for $n\lesssim 10^{3}$ in double precision, while the former algorithm was
tested for large $n$ (up to $10^{4}$ in our experiments) and demonstrated
better performance and accuracy compared to the FFT-based algorithm.
