Internal functions

recurrence_step1!(w::HCalculator)

Step 1 of the $H$ recursion: set $H^{0}_{0,0} = 1$ for every rotor, in the lower axis buffer of w, h⃗ˡ(w), whose order must be 0.

recurrence_step2!(w::HCalculator)

Step 2 of the $H$ recursion: compute the axis $H^{n}_{0,m}$, for $0 ≤ m ≤ n$, from $H^{n-1}_{0,m}$, for every rotor, where $n$ is the order of the upper axis buffer of w, h⃗ˡ⁺¹(w), and $n-1$ that of the lower one, h⃗ˡ(w). This is the recurrence of Xing et al. (2020) for the normalized associated Legendre functions, in the notation of Gumerov and Duraiswami's step 2.

recurrence_step3!(w::HCalculator{Int})

Step 3 of the $H$ recursion, for integer indices: compute the row $H^{ℓ}_{1,m}$ of the wedge of w, for $1 ≤ m ≤ ℓ$, from the axis $H^{ℓ+1}_{0,m}$ in the upper axis buffer, h⃗ˡ⁺¹(w), for every rotor. The coefficients that depend on a single integer are read from the constant table of square roots that step 2 also reads. For half-integer indices, recurrence_seed! takes its place.

recurrence_seed!(w::HCalculator{HalfOddInteger})

The half-integer form of step 3 of the $H$ recursion: compute the rows $H^{J}_{±1/2,m}$ of the wedge of w, for $1/2 ≤ m ≤ J$, from the axis $H^{j}_{0,k}$ of the integer order $j = J - 1/2$ in the lower axis buffer, h⃗ˡ(w), for every rotor. These are the two rows from which steps 4 and 5 begin; for integer indices, those rows are the row $m' = 0$, copied from the axis, and the row $m' = 1$ of step 3.

recurrence_step4!(w::HCalculator)

Step 4 of the $H$ recursion: compute the rows $H^{ℓ}_{m'+1,m}$ of the wedge of w, for $m'$ from $1$ (or $1/2$ for half-integer indices) up to $m'_{\mathrm{max}} - 1$ and $m'+1 ≤ m ≤ ℓ$, from $H^{ℓ}_{m'-1,m}$, $H^{ℓ}_{m',m-1}$, and $H^{ℓ}_{m',m+1}$, for every rotor. The first rows read are $m' = 0$ and $1$ for integer indices, and $m' = ±1/2$ for half-integer ones.

recurrence_step5!(w::HCalculator)

Step 5 of the $H$ recursion: compute the rows $H^{ℓ}_{m'-1,m}$ of the wedge of w, for $m'$ from $0$ (or $-1/2$ for half-integer indices) down to $m'_{\mathrm{min}} + 1$ and $1 - m' ≤ m ≤ ℓ$, from $H^{ℓ}_{m'+1,m}$, $H^{ℓ}_{m',m-1}$, and $H^{ℓ}_{m',m+1}$, for every rotor.

spinor_phases(R::AbstractQuaternion, [F])

Return (eⁱᵝ, z₊, z₋, cβ½, sβ½) for the rotor R, where $β$ is the Euler angle, $z₊ = e^{i(α+γ)/2}$, $z₋ = e^{i(α-γ)/2}$, $cβ½ = \cos(β/2)$ and $sβ½ = \sin(β/2)$. These are the quantities the Wigner recurrences need: $e^{i(m′α + mγ)} = z₊^{m′+m} z₋^{m′-m}$, with integer exponents even for half-integer $m′, m$, while the half-angle pair seeds the half-integer recurrence. At the poles $β ∈ \{0, π\}$ the undefined phase is set to 1; the corresponding $d$ elements vanish, so the choice is immaterial.

The half-angles are taken as $(\sqrt{W²+Z²}, \sqrt{X²+Y²})/\|R\|$, which is accurate near both poles and is non-negative, so $β ∈ [0, π]$; a rotor's double-cover sign is encoded entirely in z₊ and z₋, giving $𝔇(-R) = -𝔇(R)$ for half-integer indices.

Callers that need only the first few outputs may drop the rest: eⁱᵝ, z₊, z₋ = spinor_phases(R, F).

R need not be normalized.

These phases are not differentiable at $β = 0$ or $β = π$, where √b or √a is taken of an exact zero, and a derivative taken through them there by automatic differentiation is NaN. No local rule can repair this — sβ½ * z₋ is smooth in the rotor, but sβ½ and z₋ separately are not, so treating the zero as exact would silently drop the first-order term. Near a pole the derivatives are finite but inaccurate, the $k$-th by about $ε r^{-k}$ relative to their size at a distance $r$ from it. The values are accurate at every rotor, the poles included, so the calculators of $𝔇$ and of ${}_sY_{ℓ,m}$, which are smooth there, are never differentiated through this: the rules for automatic differentiation give their derivatives in terms of their values (see src/derivatives/kernels.jl). $H$ of a rotor has no such rules, and keeps the NaN. $d$ of a rotor is differentiated through $β$, by the same rules as $d$ of an angle. Since $β$ has a cone-shaped singularity at each pole, the derivatives of those elements of $d$ that actually have no derivative there are NaN, and those of the others are zero.

The optional second argument is the real type the phases are computed in; it defaults to float(eltype(R)). Pass the calculator's type whenever that is more precise than the rotor's own — otherwise every later step inherits the rotor type's precision.

half_angles(eⁱᵝ)

The pair $(\cos(β/2), \sin(β/2))$ for the branch $β ∈ (-π, π]$ of the phase $e^{iβ}$, computed without cancellation at either pole.

A bare phase fixes $β$ only modulo $2π$, so for half-integer indices this fixes $d$ only up to the double-cover sign $(-1)^{2ℓ}$; pass the angle $β$ itself or a Rotor if that matters.

axis_ℓ(w::HCalculator, ℓ)

Label of the integer axis that seeds the wedge of order ℓ: $ℓ$ itself for integer indices, and $ℓ - 1/2$ for half-integer ones. The axis buffers are always labelled by this Int, never by ℓ.

wedge_source(m′, m, m′ₘₐₓ)

Return (a, b, σ) such that $H_{m′,m} = σ H_{a,b}$, where (a, b) lies in the stored wedge $b ≥ |a|$, $|a| ≤ m′ₘₐₓ$.

This is the definition, for the batched engine, of how the symmetries $H_{m′,m} = H_{-m,-m′} = σ H_{m,m′} = σ H_{-m′,-m}$ of the $H$ matrix supply an element outside the stored wedge, and a read of such an element should go through it or through wedge_value. The Wigner materialize! applies the same cases a whole run of elements at a time. The harmonic one reads through this function the elements $|m| ≤ |s|$ of each spin row, or every element when it writes a block a mode at a time. Both are compared with wedge_value element by element in the tests.

An ArgumentError is thrown if no stored element can supply the requested one, which happens only when both $|m′|$ and $|m|$ exceed m′ₘₐₓ.

wedge_offset(H::HWedge, a, b)
wedge_offset(H::HWedge, a, b, m′ₘᵢₙ)

Zero-based linear offset of the first rotor's element $H_{a,b}$ in parent(H), for a stored wedge element ($b ≥ |a|$). Element iᵣ is at parent(H)[wedge_offset(H, a, b) + iᵣ].

The four-argument form takes m′ₘᵢₙ(H) from the caller, which hoists it out of a loop.

transpose_sign(m′, m)

Sign $σ$ relating the transposed element of the $H$ matrix to the original: $H_{m,m′} = σ H_{m′,m}$ and $H_{-m′,-m} = σ H_{m′,m}$. For integer indices $σ ≡ 1$; for half-integer indices $σ = sgn(m) sgn(m′)$, with $sgn(0) = 1$; see the notes on the $H$ recursion.

Which rule applies is settled by the index type, so each specialization compiles to a constant or to two cheap comparisons.

HAxis{RT}
HAxis(RT, Nᵣ, ℓₘₐₓ)

The $m'=0$, $m≥0$ axis $H^n_{0,m}$, $m = 0, …, n$, of the Hˡ matrix used in calculation of the Wigner D and d matrices, for Nᵣ rotors at once. The order $n$ is always an integer: the two axes inside an HCalculator hold successive orders, and for half-integer $ℓ$ they hold the axis of the integer order $j = ℓ - 1/2$.

• RT is the real type of the elements.

The storage parent(h) is a 1-dimensional vector holding [H[iᵣ, 0, m] for m ∈ 0:n for iᵣ ∈ 1:Nᵣ], with the rotors innermost, as in an HWedge, and the recurrence reads and writes it by that linear index, h[i]. It is allocated for orders up to ℓₘₐₓ, which maxℓ(h) reports; h.ℓ = n sets the current order, which must be an Int in 0:ℓₘₐₓ, and no other field may be changed.

This is scratch space of the recurrence inside an HCalculator, which holds two of them, and is not part of the public interface.

index_argument_error(f, names, values, integer_only=false, hint="")

The ArgumentError thrown when the index arguments of a call to f, a function defined with @index_methods, are not all Ints or all half-odd-integers. The arguments names and values are tuples of the names and the values of every index argument of the call, and each is listed in the message with its type. Below each argument that could not be an index on its own — a narrow, unsigned or wide integer, a Bool, a Rational whose integer type is not Int or whose denominator is not 2, or a range that is not a UnitRange written a:b, such as 0:2:4, 0:1:4 or Base.OneTo(4) — a sentence says what is wrong with it and how to write it instead. A call whose arguments are each acceptable but of different kinds is told which are integers and which are half-odd-integers.

With integer_only the message says that f accepts integer indices only, and a half-odd-integer argument is marked as the problem. A non-empty hint is then appended, on a line of its own, to the message of a call that has a half-odd-integer argument; it is meant to say what to use instead, such as the functions that do accept half-integer indices.

f is printed with print, so it may be the name of the function as a String, or the function or callable object itself.

δ²(ℓ, m)

$(ℓ - m)(ℓ + m + 1) = (2 d^m_ℓ)^2$, four times the square of Gumerov and Duraiswami's coefficient $d^m_ℓ = \mathrm{sgn}(m) \sqrt{(ℓ-m)(ℓ+m+1)} / 2$; the factor cancels in each recurrence, and the sign is applied separately (see sgn). Both factors are Integers for integer and half-integer indices alike, so this is exact, and the value handed to sqrt is the same one a purely integer implementation computes.

sgn(m)

Eq. (44) in Gumerov and Duraiswami (2015). Note that they define sgn differently from the usual definition — including from Julia's sign — at 0, where this is $+1$.

ϵ(m)

Eq. (7) in Gumerov and Duraiswami (2015): $ε(m) = (-1)^m$ for $m > 0$, and $1$ otherwise. The half-integer extension $ε(m) = (-1)^{⌊m⌋}$ for $m > 0$ is the unique one that keeps the $H$-form of the $m′$ ladder in Gumerov and Duraiswami's shape, as step 4 of the notes on the $H$ recursion explains. The single expression below serves both index types.