Utilities

ComplexPowers(z)

Construct an iterator to compute powers of the complex phase factor $z$, which must have magnitude approximately 1. The iterator will return the complex number $zᵐ$ for each integer $m = 0, 1, 2, \ldots$.

The parameters of the iterator's type, ComplexPowers{T, RT}, are as follows:

• T is the complex type of the powers.
• RT is the real type of their components.

A real z or one with integer components is converted to a complex floating-point number first, and a z whose magnitude is not approximately 1 is refused with a DomainError.

Example

julia> cp = ComplexPowers(cis(0.1));julia> first(cp, 5)  # Get the first 5 values from the iterator5-element Vector{ComplexF64}:                1.0 + 0.0im 0.9950041652780258 + 0.09983341664682815im 0.9800665778412417 + 0.19866933079506122im 0.9553364891256061 + 0.2955202066613396im 0.9210609940028851 + 0.3894183423086505imjulia> cis(0.1).^(0:4)5-element Vector{ComplexF64}:                1.0 + 0.0im 0.9950041652780258 + 0.09983341664682815im 0.9800665778412417 + 0.19866933079506124im 0.9553364891256062 + 0.2955202066613396im 0.9210609940028853 + 0.3894183423086506im

Notes

Stoer and Bulirsch [16] described the basic algorithm on page 24 (Example 4), though there is a dramatic improvement to be made. The basic idea is a recurrence relation, where $zᵐ$ is computed from $zᵐ⁻¹$ by adding a small increment $δz$, which itself is updated by adding a small increment given by $zᵐ$ times a constant $τ$.

Although this algorithm is numerically stable, we can improve its accuracy by factoring out $ϕ$, the power of $i$ that rotates $z$ into the sector $-π/4 < \arg z ≤ π/4$, which is as close to 1 as such a rotation can bring it. This power of $i$ can be separately exponentiated exactly and efficiently because it is exactly representable as a complex integer, while the error in the computation of $zᵐ$ is reduced significantly for certain values of $z$.

This is the recurrence of complex_powers!, which gives the same values; beyond the first few powers, it achieves a worst-case accuracy of about $0.6 m ϵ$ for $zᵐ$ — where $ϵ$ is the precision of the type of the input argument — across the range of inputs $z = e^{iθ}$ for $θ ∈ [0, 2π]$, as described under complex_powers. The original algorithm can be far worse for values of $θ$ close to $π$ — often orders of magnitude worse.

complex_powers!(zpowers, z)

Compute integer powers of z from z^0 through z^m, recursively, where m is one less than the length of the input zpowers vector.

Note that z is assumed to be normalized, with complex amplitude approximately 1; this is not checked. The algorithm, and its accuracy, are described under complex_powers.

See also: complex_powers, ComplexPowers

complex_powers(z, m)

Compute integer powers of z from z^0 through z^m, recursively, and return them as a vector of length m+1.

The number z must have complex amplitude approximately 1, or a DomainError is thrown. A real z or one with integer components is converted to a complex floating-point number first, so that the result is, for example, a Vector{ComplexF64} for z = im. The largest power m may be any non-negative integer.

This algorithm is mostly due to Stoer and Bulirsch in "Introduction to Numerical Analysis" (page 24) — with a little help from de Moivre's formula, which is essentially exp(iθ)ⁿ = exp(inθ), as well as my own alterations to deal with different behaviors in different quadrants.

There isn't usually a huge advantage to using this specialized function. If you just need a particular power, it will generally be far more efficient and just as accurate to compute either exp(iθ)ⁿ or exp(inθ) explicitly. However, if you need all powers from 0 to m, this function is several times faster than the first of those options, and about twice as fast as the second, for large m. Like those options, this function is numerically stable: measured against the exact powers of the number z it was given, its error in $zᵐ$ grows linearly in $m$, and beyond the first few powers, whose errors are a few roundings, it is at most about $0.6 m ϵ$, where $ϵ$ is the machine precision of z, and typically about a third of that, with the largest errors occurring for phases near the odd multiples of π/4. That supposes that $|z|$ is 1 to within about $ϵ/3$, as it is for a correctly rounded phase such as cis(θ); a larger departure $δ = |z| - 1$ adds up to about $m |δ|$ to the error.

The powers of -z and of im * z are exactly $(-1)ᵐ$ and $iᵐ$ times those of z, and, except where abs(real(z)) == abs(imag(z)), the powers of conj(z) are exactly the conjugates of those of z. In each case the values compare equal with ==, though the signs of zeros may differ.

See also: complex_powers!, ComplexPowers

Yindex(ℓ, m)
Yindex(ℓ, m, ℓₘᵢₙ)

Index of the mode weight $(ℓ, m)$ in the canonical ordering [f(ℓ, m) for ℓ ∈ ℓₘᵢₙ:ℓₘₐₓ for m ∈ -ℓ:ℓ], which is $ℓ(ℓ+1) - ℓₘᵢₙ^2 + m + 1$.

As for Ysize, the indices may be integers of type Int, or half-odd-integers passed as Rational{Int}s with denominator 2 or as HalfOddIntegers, and all of the indices in one call must be of the same kind. The formula gives a whole number for either kind, because for half-odd indices the quarters in $ℓ(ℓ+1)$ and $ℓₘᵢₙ^2$ cancel. ℓₘᵢₙ defaults to the smallest $ℓ$ of the given kind, which is 0 for integers and 1/2 for half-odd-integers. No bounds are checked.

See also Ysize and Yrange.

Yrange(ℓₘₐₓ)
Yrange(ℓₘᵢₙ, ℓₘₐₓ)

Vector of the $(ℓ, m)$ pairs in the canonical ordering, so that Yrange(ℓₘᵢₙ, ℓₘₐₓ)[i] is the pair stored at index i.

As for Ysize, the indices may be integers of type Int, or half-odd-integers passed as Rational{Int}s with denominator 2 or as HalfOddIntegers, and both indices in one call must be of the same kind. For half-odd-integers the pairs are of HalfOddIntegers — whatever the type of the arguments — as is every other index value the package hands back; use Rational(x) to convert an element to the more familiar Rational. ℓₘᵢₙ defaults to the smallest $ℓ$ of the given kind, which is 0 for integers and 1/2 for half-odd-integers.

See also Ysize and Yindex.

Ysize(ℓₘₐₓ)
Ysize(ℓₘᵢₙ, ℓₘₐₓ)

Total number of mode weights $(ℓ, m)$ with $ℓₘᵢₙ ≤ ℓ ≤ ℓₘₐₓ$ and $-ℓ ≤ m ≤ ℓ$, which is $(ℓₘₐₓ+1)^2 - ℓₘᵢₙ^2$.

The indices may be integers of type Int, or half-odd-integers passed either as Rational{Int}s with denominator 2 (7//2) or as HalfOddIntegers. Both indices in one call must be of the same kind; a call that mixes them, such as Ysize(0, 7//2), or that passes an integer of another type, such as an Int8, throws an ArgumentError saying so. The formula holds for either kind — for half-odd $ℓₘᵢₙ$ and $ℓₘₐₓ$ the quarters in the two squares cancel — and the result is always a whole number. ℓₘᵢₙ defaults to the smallest $ℓ$ of the given kind, which is 0 for integers and 1/2 for half-odd-integers. An ArgumentError is thrown if $ℓₘᵢₙ < 0$ or if $ℓₘₐₓ < ℓₘᵢₙ - 1$ (the empty range $ℓₘₐₓ = ℓₘᵢₙ - 1$ has size 0).

See also Yindex and Yrange.

ModeWeights(data, s=0; ℓₘᵢₙ=abs(s))
ModeWeights(data, s, ℓₘᵢₙ, ℓₘₐₓ)
ModeWeights{T}(undef, s, ℓₘᵢₙ, ℓₘₐₓ)
ModeWeights{T}(undef, s, ℓₘₐₓ)
ModeWeights(w::ModeWeights; ℓₘᵢₙ=abs(spin(w)), ℓₘₐₓ=ℓₘₐₓ(w))

Vector of mode weights $f_{ℓ,m}$ of a spin-weighted function $f = \sum_{ℓ,m} f_{ℓ,m}\, {}_sY_{ℓ,m}$, stored in the canonical ordering [f(ℓ, m) for ℓ ∈ ℓₘᵢₙ:ℓₘₐₓ for m ∈ -ℓ:ℓ] (see Yindex), together with the spin weight s and the range of $ℓ$. The parameters of the type ModeWeights{T, IT, V} are as follows:

• T is the number type of the weights.
• IT is the index type, Int or HalfOddInteger.
• V is the type of the storage, a vector of T.

A ModeWeights is not an AbstractVector; array_view gives the flat 1-based storage, which is what the transforms and the operator matrices take. It behaves as the vector of the weights in other respects: w[i], for a position i in storage order, reads or writes one weight, w[:], w[r] and w[v], for a range r or a vector v of positions, read or write several, and iteration, length, keys, sum, maximum and the other reductions see the weights one by one.

The arithmetic of mode weights as the weights of functions keeps the labels: a + b and `a

• b(and their broadcast forms) when the labels agree,complex.(a, b)` from real and

imaginary parts whose labels agree, and products and quotients with numbers, or elementwise with a plain vector of factors (a diagonal operator, such as a filter). So do zero(w), fill!(w, x), rmul!(w, a), lmul!(a, w), I * w, copyto!(w, v) from a plain vector of one weight per mode, and copyto!, copy!, axpy! and axpby! between weights whose labels agree. Anything else — adding a number to every weight, the elementwise product of two sets of weights, conj.(w), abs2.(w) — would label numbers that are not the mode weights of any function, and is an error; arithmetic on the raw numbers goes through array_view(w). A broadcast whose result is labelled is a ModeWeights, and one that would extend the weights to a longer vector is refused. map returns plain numbers, and so do similar(w, n) and similar(w, T, dims), which give storage of the size asked for, while similar(w) and similar(w, T) keep the labels.

dot(a, b) is the inner product of the two functions, and refuses two ModeWeights whose labels differ. w' and transpose(w) are those of the storage, so a' * b is a product with a plain matrix, which is refused, as described below. a == b, isequal(a, b) and a ≈ b compare the labels as well as the numbers, and are false when the labels differ; against a plain vector they compare the numbers alone, and so hash(w) is the hash of the numbers. In addition

• w[ℓ, m] reads or writes the weight of mode $(ℓ, m)$,
• w[ℓ, :] is a DegreeBlock view of the weights for one $ℓ$, indexed by m ∈ -ℓ:ℓ,
• modes(w) is the vector of (ℓ, m) pairs in storage order,
• spin(w), ℓₘᵢₙ(w), ℓₘₐₓ(w) are the parameters, and parent(w) is the storage,
• the differential operators L², Lz, L₊, L₋, Lx, Ly, R², Rz, R₊, R₋, ð, ð̄ give a new ModeWeights with the spin weight adjusted where appropriate and the range of $ℓ$ of w, written either ð * w or ð(w). These build no matrix: the operator is applied by a loop, so the only allocation is the result, and mul!(w′, ð, w) into a correctly labelled destination allocates nothing at all (w′ may also be a bare vector at least as long as the result, which then comes back as a ModeWeights over it),
• a product with a plain matrix, A * w, is refused, because nothing about a matrix of numbers says what spin weight or range of $ℓ$ it was built for, so the labels could not be checked: A * array_view(w) is the product with the raw numbers, sYlm(R⃗, ℓₘₐₓ, s) * w evaluates the function at rotors, and op * w or mul!(w′, op, w) applies a differential operator with the labels checked, and
• w(R) evaluates the function at the rotor R (see sYlm).

When constructed from data alone, ℓₘₐₓ is deduced from length(data) and ℓₘᵢₙ, which must match exactly. The data vector is used as storage, not copied, so it must keep its length for as long as the ModeWeights is in use; the operators, w[ℓ, :], and w[ℓ, m] outside @inbounds throw a DimensionMismatch when they find that it has been resized. The undef forms allocate uninitialized storage of type T instead, and in the three-argument form ℓₘᵢₙ defaults to abs(s), as it does when data is given. The keyword ℓₘᵢₙ may also be spelled ell_min. show writes the labels, as ModeWeights{ComplexF64} with s=2, ℓ ∈ 2:8, and the text/plain form adds the weights.

Changing the range of $ℓ$

The spin-changing operators keep the range of $ℓ$ of their input, so the result of ð * w for w with $ℓₘᵢₙ = |s|$ starts at $ℓ = |s|$, not at $|s + 1|$ as the weights of its new spin weight usually do. Weights of one spin weight but different ranges of $ℓ$ cannot be added or compared, and the destination of mul! must have exactly the labels of the result. ModeWeights(w; ℓₘᵢₙ, ℓₘₐₓ) copies w into another range of $ℓ$, with the same spin weight: it fills with zeros every $ℓ$ that w does not hold, drops the entries with $ℓ < |s|$, which belong to no harmonic, and drops those above the ℓₘₐₓ requested. With the defaults, it copies w into $|s| ≤ ℓ ≤ ℓₘₐₓ(w)$, the range the transforms' analysis produces, so it is a copy only when w already has that range. The keywords may also be spelled ell_min and ell_max, and must be indices of the kind of w's own.

w₂ = ModeWeights(ð * (SSHT(1, ℓₘₐₓ) \ f); ℓₘᵢₙ=2)   # the spin-2 weights on ℓ ∈ 2:ℓₘₐₓw₂ ≈ SSHT(2, ℓₘₐₓ) \ g

Several functions at once

A transform applied to several functions at once, as the columns of a matrix, returns a plain array whose columns are the mode weights, and a ModeWeights holds one vector of them. ModeWeights(view(W, :, j), s) labels column j without copying it, so that an operator or a rotation can be applied column by column with no allocation:

for j ∈ axes(W, 2)    mul!(view(W′, :, j), ð, ModeWeights(view(W, :, j), s))end

The weights are real or complex numbers, and a quaternion element type is refused. A Rotor is a quaternion, so R * w is an error; the rotation of the function by R is D(R, ℓₘₐₓ(w)) * w (see D).

Half-integer indices

The spin weight and the range of $ℓ$ may be half-integers, passed as Rationals with denominator 2 — as in ModeWeights(data, 1//2) or ModeWeights{T}(undef, 1//2, 1//2, 7//2) — or as HalfOddIntegers, in which case every $ℓ$ and $m$ of the ordering is a half-odd-integer, and ℓₘᵢₙ may be as small as 1//2. The indices in one call must all be of one kind, integers of type Int or half-odd-integers; a call that mixes them, such as ModeWeights(data, 1//2, 0, 7//2), is an error, as is an integer index of another type, such as an Int32, which is to be converted with Int first. The parameters are stored as HalfOddIntegers, which is also what modes(w) and the axis of w[ℓ, :] are made of; w[ℓ, m] accepts either spelling. For such a w, w[ℓ, :] is a DegreeBlock, indexed by m ∈ -ℓ:ℓ, exactly as it is for integer indices. The natural indices of w[ℓ, m] and w[ℓ, :] obey the same rules as those of the constructors: an integer w takes Ints, and an index of another kind or integer type is refused with the reason.

ModeWeights(w::ModeWeights; ℓₘᵢₙ=abs(spin(w)), ℓₘₐₓ=ℓₘₐₓ(w))

A copy of the mode weights w over the range ℓₘᵢₙ:ℓₘₐₓ of $ℓ$, with the same spin weight, zero where w holds no weight.

Every $ℓ$ of the new range that w does not hold is filled with zeros, the entries of w with $ℓ < |s|$, which belong to no harmonic, are dropped, and so are those above the ℓₘₐₓ requested. The defaults give the range $|s| ≤ ℓ ≤ ℓₘₐₓ(w)$, which is what the transforms' analysis produces. The result has new storage, of the element type of w, even when the range is that of w itself. The keywords may also be spelled ell_min and ell_max; they must be indices of the kind of w's own, integers of type Int or half-odd-integers, each a HalfOddInteger or a Rational{Int} with denominator 2.

This is the way to bring weights of one spin weight but different ranges of $ℓ$ together: the spin-changing operators keep the range of their input, while +, ≈, dot and the destination of mul! require the labels to agree exactly.

w = SSHT(1, 8) \ f                  # s = 1, ℓ ∈ 1:8ðw = ModeWeights(ð * w; ℓₘᵢₙ=2)      # s = 2, ℓ ∈ 2:8, as SSHT(2, 8) \ g isð̄w = ModeWeights(ð̄ * w; ℓₘᵢₙ=0)      # s = 0, ℓ ∈ 0:8, with ℓ = 0 filled with zeros

w[ℓ, :]

A view of the mode weights of the ModeWeights w for the given $ℓ$, indexed by m ∈ -ℓ:ℓ. This is a DegreeBlock for either kind of index; where the indices are half-odd-integers they may be passed either as Rationals or as HalfOddIntegers. Writing through the view writes into w.

w[ℓ, m]

The mode weight of $(ℓ, m)$ in the ModeWeights w. For a w with half-integer indices, ℓ and m may be passed as Rational{Int}s with denominator 2 — w[3//2, 1//2] — or as HalfOddIntegers, either way for each; for a w with integer indices they must be Ints. An index of another kind or type is refused with an ArgumentError that says why, as the constructors refuse one.

modes(w::ModeWeights)

The (ℓ, m) pairs of w, in storage order (see Yrange).

check_rotor_type(calc, data)

Throw an ArgumentError unless data would give the element type calc already works in.

A calculator's element type is fixed by the data it was constructed from, so replacing that data later — through set_R!, set_β!, set_θ! or similar(calc, data) — cannot change it. Rather than convert silently, which is how precision gets lost without anyone choosing to lose it, the mismatch is an error and the caller converts whichever side they meant. Both sides are given by floattype, whose method for each calculator is defined alongside it.

not_a_rotor(R)

The message for rotor data R, or for rotor data of type R, that is a quaternion but does not denote a rotation, which is to say a QuatVec, or a vector of them.

These functions are defined on the rotation group, so a rotation is what they take: a Rotor, or any other Quaternion, which denotes the rotation of its normalization, since the recurrence divides out its magnitude. A QuatVec is a different thing — it represents a vector, and reading one as a rotation by $π$ about its own direction would be a category error rather than a convenience; exp(v/2) gives the rotation that a vector generates.

nrotors(R)

The number of rotors Nᵣ that the rotor data R describes: one for a single rotor, angle or phase, and length(R) for a vector of them. This is how a calculator learns its batch size, which is why there is no Nᵣ keyword argument on any constructor.