Module sh 

Real spherical harmonics, evaluated directly from Cartesian components.

Conventions

Orthonormal real spherical harmonics. The CS const parameter selects the phase convention: false for the standard real-SH tables (sphericart, most math references), true for the Condon-Shortley phase. Everything below describes CS = false. See the phase section for what the other one changes.

\int_{S^2} Y_{\ell m}^2 \, d\Omega = 1,
\qquad
Y_{\ell m} =
\begin{cases}
\sqrt{2}\, K_\ell^m P_\ell^m(\cos\theta)\cos(m\varphi) & m > 0 \\
K_\ell^0 P_\ell(\cos\theta) & m = 0 \\
\sqrt{2}\, K_\ell^{|m|} P_\ell^{|m|}(\cos\theta)\sin(|m|\varphi) & m < 0
\end{cases}

with $K_\ell^m = \sqrt{\tfrac{2\ell+1}{4\pi}\tfrac{(\ell-m)!}{(\ell+m)!}}$ and $P_\ell^m$ the associated Legendre functions without $(-1)^m$. So $Y_{00} = \sqrt{1/4\pi}$, $Y_{1,-1} = \sqrt{3/4\pi}\,y$, $Y_{10} = \sqrt{3/4\pi}\,z$, $Y_{11} = \sqrt{3/4\pi}\,x$.

Outputs are written in the flat l * (l + 1) + m order (m from -l to l), the layout every SH-lighting pipeline uses.

The Condon-Shortley phase

The two conventions differ by $(-1)^{|m|}$: odd |m| is negated, even |m| agrees exactly. Which one a body of data was projected against is not recoverable from the data (the difference is invisible in any rotationally-averaged or squared quantity), so mixing them is a silent, plausible-looking wrong answer. Hence the explicit parameter rather than a fixed choice.

Sloan’s widely-copied SHEval generated code (Efficient Spherical Harmonic Evaluation, JCGT 2(2), 2013) does carry the phase. Its diagonal recurrence is P_m^m = (1 - 2m) P_{m-1}^{m-1}, negative for every m >= 1. Its order-3 listing emits pSH[3] = -0.48860251 * x, matching CS = true here, while CS = false gives +0.48860251 * x.

CS is baked into the constant table, so neither choice costs an instruction. Internally it is applied in two places, the second easy to overlook: the diagonal seeds for odd m (which propagates to a whole column, and to both the +m and -m slots that share it), and every z-derivative ratio in f, because that ratio crosses between adjacent columns whose signs always disagree.

Algorithm

No trigonometry and no division anywhere. The evaluation factors each harmonic as $Y_{\ell,\pm m} = q_\ell^m(z) \cdot \{c_m, s_m\}$ where

• $c_m + i s_m = (x + iy)^m$, accumulated by the complex-multiplication pair recurrence. Since $x + iy = \sin\theta\, e^{i\varphi}$ on the unit sphere, this is $\sin^m\theta \{\cos, \sin\}(m\varphi)$, i.e. the $\sin^m\theta$ factor of $P_\ell^m$ moved into the azimuthal term, which removes the $1/\sin\theta$ pole from every recurrence (the factoring used by sphericart, Bigi et al., J. Chem. Phys. 159, 064802, 2023).
• $q_\ell^m(z)$ is the fully-normalized sin-factored associated Legendre part, via the standard normalized three-term recurrences (Holmes & Featherstone 2002, J. Geodesy 76): a constant diagonal, one $\sqrt{2m+3}\, z$ step, then $q_\ell^m = a_\ell^m z\, q_{\ell-1}^m - b_\ell^m q_{\ell-2}^m$. All coefficients are precomputed at compile time (ShConsts). Intermediate values stay O(1), so there is no overflow at any order either format can index.

Cost is O(L^2) FMAs per call (two per harmonic past the seeds) with zero transcendentals.

Domain and gradient semantics

(x, y, z) is assumed to be a unit vector. Nothing renormalizes. Off the unit sphere the recurrences still evaluate a perfectly good polynomial in (x, y, z) (the one that agrees with $Y_{\ell m}$ on the sphere), which is exactly what sh_d_impl’s derivatives differentiate: the ambient Cartesian gradient of that polynomial form, evaluated at the given point. This is the convention machine learning interatomic potentials and finite-difference checks want. A caller who needs the tangential (spherical) gradient projects out the radial component: g_tan = g - (g . n) n.

The derivative combinations are exact identities on the recurrence outputs: $\partial_x c_m = m c_{m-1}$, $\partial_y c_m = -m s_{m-1}$ (and the mirrored pair for $s_m$), and $\partial_z q_\ell^m = f_\ell^m q_\ell^{m+1}$ where $f_\ell^m$ is a tabulated norm ratio. So the gradient pass reuses every value the value pass produced and adds no new recurrences.

Structs

ShTable

Precomputed recurrence coefficients for all (l, m) with l <= L.

Constants

MAX_DEGREE

Highest degree the stamped ladders below cover. Beyond it the kernels fall back to the rolled, runtime-coefficient path (sh_table_impl + sh_eval_impl), which is correct at any degree but roughly an order of magnitude slower. Extending the ladder is mechanical: append literals to every 0 1 2 ... 16 list.

Traits

ShConsts

Compile-time spherical-harmonic coefficient tables for one element type, at one degree and one phase convention.

Functions

sh_d_impl

sh_impl plus the ambient Cartesian gradient of every harmonic.

sh_eval_d_impl

sh_eval_impl plus the ambient Cartesian gradients. The rolled counterpart of sh_d_impl, with the same two-pass structure and gradient semantics.

sh_eval_impl

Evaluates all harmonics through degree L from a table filled by sh_table_impl.

sh_eval_lifted_impl

sh_eval_impl over a table stored in W’s primal type, each coefficient lifted through from_primal as it is read.

sh_eval_mixed_impl

sh_eval_impl with the coefficients kept in a different, simpler type than the values.

sh_impl

All real spherical harmonics through degree L at the unit direction (x, y, z).

sh_table_impl

Computes the recurrence coefficients for degree L into a runtime table.

tri

Triangular index base: q_l^m lives at tri(l) + m.