Half-integer indices

Wigner's matrices and the spin-weighted spherical harmonics are defined for half-integer $ℓ, m', m$ and $s$ as readily as for integer ones. These are the irreducible representations of $𝐒𝐩𝐢𝐧(3) ≅ 𝐒𝐔(2)$ that do not descend to $𝐒𝐎(3)$, and they are supported throughout this package — by $𝔇$, $d$ and ${}_{s}Y_{ℓ,m}$, by the mode weights and the operators acting on them, and by the transforms — at the same accuracy and essentially the same speed as integer indices.

Because half-integers appear in every part of the interface, the way they are input, the types they produce, and the few places where they behave differently are all collected on this page.

Creating a half-integer

A half-integer index can be input as a Rational{Int} whose denominator is exactly 2:

𝔇 = D(R, 7//2)𝔡 = d(β, 7//2)calculator = DCalculator(R⃗, 7//2)Y = sYlm(R, 7//2, -3//2)𝒯 = SSHT(1//2, 7//2)

Note that Rational numbers must be written with the double-slash // operator, not a single slash; 7//2 is a Rational, but Julia immediately converts 7/2 to 3.5 — a Float64 — which will produce an error. Optional keyword limits such as m′ₘₐₓ, m′ₘᵢₙ, mₘₐₓ and mₘᵢₙ (or their ASCII spellings mp_max, mp_min, m_max and m_min) may also occur in some of these functions, and follow the same rules. A HalfOddInteger, described below, is accepted wherever such a Rational is.

The indices within a single call must all be of one kind. An integer $ℓ$ goes with an integer $m$ and an integer spin weight; a half-odd $ℓ$ with half-odd ones. There is no mixed case, because no single representation contains both. A call that mixes them is refused with a specific, explanatory error message — rather than with the bare MethodError that dispatch alone would produce:

julia> sYlm(R, 7//2, 1)ERROR: ArgumentError: The indices of one call to `sYlm` must all be integers of type `Int`, like 3, or all be half-odd-integers, each a `HalfOddInteger` or a `Rational{Int}` with denominator 2, like 7//2.  This call has    ℓₘₐₓ = 7//2::Rational{Int64}    s = 1::Int64and so mixes integers (s) with half-odd-integers (ℓₘₐₓ).

The same kind of message explains every other refusal of an index, as described under Functions that take indices below.

The `HalfOddInteger` type

These indices are not kept as Rationals internally. At the boundary of every public function they are converted to HalfOddInteger, a type that stores the odd numerator of $n//2$ and nothing else. It exists for two reasons.

It makes the kind of index a property of the type. Every container and calculator is parameterized by its index type, so the integer and half-integer cases dispatch to separate, separately compiled code, rather than branching on a denominator at run time.
It keeps the recurrences in integer arithmetic. The sum or difference of two HalfOddIntegers is an Int, and so is $2ℓ$, so a recurrence coefficient such as $(ℓ-m)(ℓ+m+1)$ is computed entirely in integers while being written exactly as the references write it. The same expression in Rational arithmetic is roughly thirty times slower.

However, this means that a HalfOddInteger is not a number in the usual sense — or a Number in the specific Julia sense. As a result, the type is deliberately spare. It defines only the operations the recurrences and the public interface actually need, and in particular defines no promote_rule, so that an unanticipated operation is a MethodError rather than a silent — and invisibly slow — fall-back to Rational. Because a half-odd-integer is never zero and never one, zero, one, and oneunit throw rather than returning a value of some other type, and for the same reason a product with an odd integer, which would not be a half-odd-integer, is refused with a DomainError, while 2ℓ is simply an Int. A HalfOddInteger converts exactly to Float16, Float32, Float64 and BigFloat, as Float64(ℓ), and to the types of DoubleFloats once that package is loaded; another floating-point type T is reached through T(Rational(ℓ)).

Together with the Integers, these make up IntegerHalf, the union of index types over which the recurrences are defined. Neither name is exported, and neither needs to be constructed by hand: because the conversion happens at the interface boundary, 7//2 remains a natural thing to write everywhere — though it is also acceptable to construct a HalfOddInteger directly and pass that. What does deserve noticing is that the values coming out — from ℓ, m′ₘₐₓ, axes, and the other accessors — are HalfOddIntegers, though they display simply as 7//2 and compare equal to it. Call Rational(i) to convert one back.

What the calls return

A block with integer indices could be an OffsetArray, indexed by its natural $m$ and $m'$, but half-odd indices cannot be handled that way, because OffsetArrays cannot use non-integer axes. This package has its own family of containers instead, and they are what every call returns, whichever kind of index is in play, so there is no table of correspondences to learn:

call	what comes back
`D(R, ℓₘₐₓ)`, `d(β, ℓₘₐₓ)`	a `WignerSeries` of blocks
a block of `D`	a `WignerMatrix` of complex numbers
a block of `d`	a `WignerMatrix` of real numbers
a block of a `DCalculator` or `dCalculator` built from a vector of rotors	`WignerMatrixBatch`
`sYlm(R, ℓₘₐₓ, s)`	`HarmonicValues`
a block of an `sYlmCalculator`, one spin weight	`DegreeBlock`
the same, over a batch of rotors	`DegreeBlockBatch`
a block of a range of spin weights	`SpinMatrix`
the same, over a batch of rotors	`SpinMatrixBatch`

The reason for keeping OffsetArrays off the integer path is not just uniformity, but also safety. As is set out under Containers: an OffsetArray with non-trivial offsets accepts * and mul! and returns silently wrong answers. array_view is the explicit route to a plain 1-based array, for either kind of index, and relabel is the way back.

For half-integers, the index $ℓ$ runs over 1//2, 3//2, …, ℓₘₐₓ, so ℓₘᵢₙ is 1//2 rather than 0, and that is where a loop over a calculator begins. Indexing a series with something that is not one of those values — 𝔇[3], 𝔇[1//1] — throws an error.

These containers support [m′, m] (or [iᵣ, m′, m] where relevant), checkbounds, axes, size, size(w, d), length, ndims, eltype, parent, copy, similar, collect, Array, Matrix, ==, isequal, hash, ≈, iteration and show, and a WignerMatrixBatch additionally gives w[iᵣ] — a WignerMatrix view of one rotor's block. The comparisons count the labels as well as the numbers: two blocks are equal only if they are the same kind of block, for the same $ℓ$ and the same ranges of indices, holding equal values. The blocks are deliberately not AbstractMatrixes, so linear algebra does not apply to them directly: operations like B', B * C, lu(B) and the like are MethodErrors. Broadcasting (B .+ 1) returns an ordinary 1-based Array, dropping the natural indices, but it does check them first: every container in one broadcast, including the destination of B .= C, must have the same labels, and a mismatch is refused with an ArgumentError. Call array_view for a 1-based view of the same storage, on which BLAS works at full speed, or Matrix(B) for an independent copy.

$β$, the double cover, and the $H$ wedge

Half-integer $d$ has period $4π$ in $β$ rather than $2π$: $d^{(ℓ)}(β + 2π) = -d^{(ℓ)}(β)$. An angle or a Rotor therefore pins the result down unambiguously, but a bare phase $e^{iβ}$ fixes $β$ only modulo $2π$; in that case the branch $β ∈ (-π, π]$ is used, so d(cis(β), ℓₘₐₓ) agrees with d(β, ℓₘₐₓ) to within rounding inside that branch and differs by the sign $(-1)^{2ℓ}$ outside it.

The double cover shows up in $𝔇$ as $𝔇^{(ℓ)}(-R) = -𝔇^{(ℓ)}(R)$, which holds exactly — the two results are bitwise negatives of each other, with no rounding. A Rotor passed to d contributes only $β ∈ [0, π]$, so d(-R, ℓₘₐₓ) == d(R, ℓₘₐₓ): the sign lives entirely in the $α$ and $γ$ phases.

One caveat applies to the $H$ wedge underneath. For half-integer indices $H$ is symmetric only up to the sign $σ = \mathrm{sgn}(m)\,\mathrm{sgn}(m')$, so code that reads calc.Hˡ and transposes it by hand will get the wrong sign for half of the elements; use wedge_value, which applies $σ$ for you, and accepts the indices as HalfOddIntegers or as Rational{Int}s with denominator 2. See the notes on the $H$ recursion.

Half-integer spin weight

Two things differ mathematically when the spin weight itself is half-odd. The prefactor $(-1)^s$ in the definition of ${}_{s}Y_{ℓ,m}$ is then $\pm i$; this package uses the principal branch $(-1)^s ≡ e^{iπs} = i^{2s}$, settled in the conventions summary. Consequently ${}_{s}λ_{ℓ,m}(θ)$ — the value at $(θ, ϕ=0)$, which is real for integer $s$ — is not real here.

The other difference is in what a function of half-integer spin weight is a function of. The two rotors $±𝐑$ above a point of the sphere give values of opposite sign, so "the value at $(θ, ϕ)$" is defined only once the rotor is named. That has consequences for the sampling that the transforms do, and is discussed on the transformations page.

What is not supported

Half-integer indices are available throughout most of the package, except where otherwise noted. Three things do not accept them. Ylm and YlmCalculator are integer-only, because the ordinary spherical harmonics are the spin-weight-zero members of the family, and no half-integer series contains $s = 0$; a half-integer $ℓₘₐₓ$ is refused with a message pointing to sYlm and sYlmCalculator. The two equiangular grids, driscoll_healy_pixels and mcewen_wiaux_pixels (with their _rotors counterparts), take only integer indices, both the spin weight (which they ignore) and the band limit, since no transform uses them by default; a half-integer index is refused with a message pointing to the golden-ratio, Leja, and sorted-ring pixelizations, which do accept one. And the "Minimal" transform method is integer-only, as is minimal_rings, because its bookkeeping of rings and aliased modes is written for integer indices and has not (yet) been extended; SSHT(1//2, 7//2; method="Minimal") is refused with a message naming the "RS" and "Matrix" methods, which do accept half-integer indices.

Functions that take indices

The public functions and constructors that take indices are defined with @index_methods, which is what applies the rules on this page and writes the messages of the refusals. The index arguments of such a definition are annotated with one of three markers: IndexType for a single index, IndexRange for a range of them, and IndexOrRange for either, as the spin weight of sYlm may be. From that one definition the macro generates the methods that convert a Rational{Int} with denominator 2 to a HalfOddInteger, that refuse every other value with an explanation, and that pass on to the body of the definition indices which are all Ints or all HalfOddIntegers.

The integer indices are Ints. A narrow integer such as an Int8, an unsigned integer, a Bool, and a wide one such as a BigInt or an Int128 are all refused, with a message saying how to write the index instead, because the index arithmetic is not closed under those types: $ℓ^2$ overflows an Int8 at $ℓ = 12$, and $-m$ wraps around in a UInt. Likewise, a half-odd-integer given as a Rational must be a Rational{Int}, a Rational that is a whole number, such as 2//1, must be written as the integer, and any other Rational, such as 5//3, is not an index at all. A range of indices must be a UnitRange, written a:b, so that 0:2:4 is refused where 0:4 is accepted; so are 0:1:4 and Base.OneTo(4), which have the right elements but are of other types, with a message saying to write them as 0:4 and 1:4. The markers and the macro are public, so that a package building on this one can define its own functions of indices in the same way, and have them accept exactly what these accept.

Two kinds of call are written by hand instead, because the kind of index they accept is fixed by an object that has already been built, rather than by the other arguments of the call. Indexing a container, as in w[ℓ, m], 𝔇[ℓ] or 𝔇[ℓ][m′, m], follows the rules above, and in addition requires each index to be of the container's own kind; a position in the flat storage, as in w[i], is an ordinary array index. The other kind, recurrence!(calc, ℓ) and wedge_value(Hˡ, iᵣ, m′, m), follows the same rules, and requires its indices to be of the kind that the calculator or the wedge was built for.

Docstrings

SphericalFunctions.HalfOddInteger — Type

HalfOddInteger(x)

A half-odd-integer — one of $…, -3//2, -1//2, 1//2, 3//2, …$ — stored as its odd numerator, so that the represented value is numerator//2. Together with the Integers these make up the IntegerHalfs, the index types over which the Wigner recurrences are defined.

The point of the type is that the arithmetic which actually occurs in those recurrences lands back in Int:

julia> using SphericalFunctions: HalfOddIntegerjulia> ℓ = HalfOddInteger(7//2)7//2julia> 2ℓ + 1  # an `Int`; compiles to `ℓ.numerator + 1`8julia> m = HalfOddInteger(3//2);julia> (ℓ - m) * (ℓ + m + 1)  # an `Int`12

Sums and differences of two HalfOddIntegers are Ints; adding or subtracting an integer gives back a HalfOddInteger. Multiplication is defined only by an even integer, and returns an Int; 2ℓ is the case that matters and it costs nothing. Multiplying by an odd integer throws a DomainError. An integer operand of any type is converted to Int first, so that one too large for an Int throws an InexactError rather than wrapping around.

A HalfOddInteger is ordered against integers and against other HalfOddIntegers, and is == to the Rational or the floating-point number of the same value, with which it also hashes alike. It is never == to an integer, and isinteger is false for it.

The type is deliberately spare: it defines only the operations the recurrences and the public interface actually need, and in particular it defines no promote_rule. An operation that has not been anticipated is therefore a loud MethodError rather than a silent fall-back to Rational arithmetic, which is roughly thirty times slower and would otherwise be invisible. If you need an operation that is missing, add it deliberately.

Because a HalfOddInteger is never zero and never one, zero, one and oneunit all throw rather than returning a value of some other type.

Users need not construct these directly: every public entry point in this package that accepts a half-integer index also accepts a Rational{Int} with denominator 2, and converts it automatically. Use Rational(x) to convert back, or Float64(x) for a floating-point value. The conversion to a float type is defined for Float16, Float32, Float64 and BigFloat, and for Double16, Double32 and Double64 when DoubleFloats is loaded; any other float type T converts through T(Rational(x)).

See also IndexType and IndexRange.

source

SphericalFunctions.@index_methods — Macro

@index_methods function f(args...; kwargs...) ... end
@index_methods integer_only function f(args...; kwargs...) ... end
@index_methods integer_only "hint" function f(args...; kwargs...) ... end

Define f for index arguments of any type a caller may reasonably write, while the body sees only indices of type Int or only HalfOddIntegers.

The index arguments of the definition are annotated with one of the markers IndexType, IndexRange and IndexOrRange, either directly, as in ℓₘₐₓ::IndexType, or through a type variable bounded by a marker, as in ℓₘₐₓ::IT, s::IT with where {IT<:IndexType}, in which case the body may use IT. Other arguments, positional defaults and where clauses are written as usual. The macro generates four methods of f:

A fallback, whose index arguments are typed with their markers, and which throws the ArgumentError of index_argument_error, naming every index argument with its value and type and saying what is wrong. This is the method that a call reaches if its indices are not all of one kind, or if any of them is an integer of a type other than Int, 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: one whose step is not 1, or one such as 0:1:4 or Base.OneTo(4), whose message says to write it 0:4 or 1:4. The docstring written above the macro call is attached to this method.
A conversion method, whose index arguments are typed Union{Rational, HalfOddInteger}, or the unit ranges of those, which converts each Rational{Int} with denominator 2 to a HalfOddInteger, and a range endpoint by endpoint, and calls f again with the keyword arguments forwarded untouched. Any other Rational is refused with the fallback's error.
A work method whose index arguments are Ints or UnitRange{Int}s (with IT<:Int for a type variable), whose keyword arguments are exactly those written, and whose body is the body written.
The same with HalfOddInteger and UnitRange{HalfOddInteger}.

Defaults, positional and keyword, are evaluated in the work methods, where the positional indices have already been converted, so that a default computed from an index, such as Nϕ::Int=2ℓₘₐₓ+1, sees an Int or a HalfOddInteger, and a default may use IT. The exception is the default of an argument that no index argument precedes, including that of the first index argument itself: no index is known where it is evaluated, so it may not use IT, and a call that omits it is completed by one more generated method, which evaluates it and calls f again. A definition with optional positional arguments has a fallback and a conversion method for each number of arguments with which it may be called, and the docstring is attached to the fallback that takes them all. A call that is refused is refused before any default is evaluated.

Keyword arguments take no part in dispatch, so each keyword annotated with a marker, or with an index type variable, which then stands for its marker, or with the Union of either with Nothing, is normalized at the start of each work method, after the defaults have been evaluated. A keyword of the kind of the positional indices is used as it is, a Rational{Int} with denominator 2 is converted where that kind is HalfOddInteger, nothing is passed through where the annotation allows it, and anything else is an ArgumentError naming the keyword. A default that refers to an earlier keyword index sees that index normalized in the same way. This is what makes the ASCII-alias pattern

@index_methods function f(ℓ::IndexType; mp_max::IndexType=ℓ, m′ₘₐₓ::IndexType=mp_max)    # only m′ₘₐₓ is used hereend

accept either spelling, of either kind. A keyword that is not annotated with a marker is passed through untouched.

With integer_only, the conversion method and the fourth method both throw an ArgumentError saying that f accepts integer indices only. The fourth method is kept, rather than omitted, so that the methods of f remain free of ambiguities when another definition of f has less specific non-index arguments. A string written after integer_only is a hint, which is appended to that message when the call has a half-odd-integer index, to say what to use instead:

@index_methods integer_only (    "Functions of half-integer spin may be sampled with `golden_ratio_spiral_pixels`, "    * "`leja_pixels` or `sorted_ring_pixels`.") function driscoll_healy_pixels(s::IndexType, ℓₘₐₓ::IndexType, ::Type{T}=Float64) where {T}    ...end

The hint is a string literal, or a concatenation of string literals with * as here, so that a long one may be written over several lines; it is fixed when the definition is expanded.

The definition must have at least one positional index argument, since that is what the work methods dispatch on, and every index argument must be named. A type variable bounded by a marker may annotate index arguments only, and one bounded by IndexOrRange may annotate only one of them, because the arguments it annotated would then all have to be of one type, all indices or all ranges. A short-form definition with both a return type and a where is written with the signature in parentheses, as (f(ℓ::IT)::IT) where {IT<:IndexType} = ..., since otherwise the where applies to the return type. Macros applied to the definition, such as @inline but not @generated, are written between @index_methods and function, and are applied to each generated method.

Example

"""    Ysize(ℓₘᵢₙ, ℓₘₐₓ)..."""@index_methods function Ysize(ℓₘᵢₙ::IT, ℓₘₐₓ::IT) where {IT<:IndexType}    ...end

source