Difficulty defining promotion rules


I’m trying to define promotion rules for a type representing multicomplex numbers (Multicomplex number - Wikipedia). These have a series of complex numbers, im1, im2, etc. Internally, these are represented as static arrays of size 2^N:

using StaticArrays
struct Multicomplex{T,N,C} <: Number
im1 = Multicomplex{Bool,1,2}(SVector{2}(false, true))
im2 = Multicomplex{Bool,2,4}(SVector{4}(false, false, true, false))

I’m trying to define a promotion rule so that Multicomplex numbers (e.g. im1) get promoted to larger sizes when needed in calculations:

             ::Type{<:Multicomplex{S,N2,C2}}) where
             {T<:Real,S<:Real,N1,C1,N2,C2} =
# constructors are also defined for conversion, e.g.
#Multicomplex{T,2,4}(m::Multicomplex{S,1,2}) where {T,S} =
#    Multicomplex{2}(SVector{4,T}(m.value[1], m.value[2], zero(T), zero(T)))

However, this rule doesn’t seem to work, as type parameters are just dropped when they don’t match in promote_type:

julia> promote_type(Multicomplex{Float64,1,2}, Multicomplex{Float64,1,2})
Multicomplex{Float64, 1, 2}

julia> promote_type(Multicomplex{Float64,1,2}, Multicomplex{Float64,1,4})
Multicomplex{Float64, 1}

julia> promote_type(Multicomplex{Float64,1,2}, Multicomplex{Float64,2,4})

julia> promote_type(Multicomplex{Float64,1,2}, Multicomplex{Bool,2,4})

Any advice would be appreciated. Thank you!

There are two problems there.

First, you need to import Base.promote_rule (your definition is being ignored by promote_type). Second, maximum is for arrays, you want max:

julia> function Base.promote_rule(::Type{<:Multicomplex{T,N1,C1}},
                    ::Type{<:Multicomplex{S,N2,C2}}) where

julia> promote_type(Multicomplex{Float64,1,2}, Multicomplex{Float64,1,4})
Multicomplex{Float64, 1, 4}

Two stupid mistakes :grinning: I always get max and maximum confused, and had imported most other functions from Base but forgot this one! Thank you for the quick response!

What does the second type parameter, N, encode? It seems redundant if it’s just C = 2^N. Can’t you use only the C parameter for simplicity?

Another thing, did you consider using Tuple instead of SVector? Or do you need StaticArray functionality?

Yes, it is encoding C=2^N. My thinking is this: I’d originally tried to avoid multiple parameters by having type {T,N} only, but then I couldn’t correctly set the size of the SVector. I suppose it’s possible to have the type as {T,C}, but most users will be thinking in terms of N and not C, so I’d rather expose this. By having the type as {T,N,C} (and constructors specified in terms of T and N), users can ignore the last parameter C and just work with Multicomplex{T,N} objects.

If you’ve a better suggestion though it would be welcome - I’m a relative newbie.

My thinking for SVector was that I’ll be doing arithmetic on the contents, e.g.

Base.:(+)(a::Multicomplex{T,N}, b::Multicomplex{S,N}) where {T,S,N} =
    Multicomplex{N}(a.value + b.value)
Base.:(-)(a::Multicomplex{T,N}, b::Multicomplex{S,N}) where {T,S,N} =
    Multicomplex{N}(a.value - b.value)

This wouldn’t be so straightforward with Tuple?

I see. It’s a bit hard to say, but it seemed like values would be constructed as

im2 = Multicomplex{Bool,2,4}(SVector{4}(false, false, true, false))

etc, so then the user would anyway have to consider the number of values to input, and think ‘in terms of C’. I guess it depends on how this will be used in the end.

Tuples can also do arithmetic with the help of broadcasting

julia> (1,2,3) .+ (4,5,6)
(5, 7, 9)

which should be as performant as staticarrays (StaticArrays.jl is built on tuples.)

One problem, at least with the code you shared, is that nothing is stopping me from writing

im1 = Multicomplex{Bool,11,2}(SVector{2}(false, true))
im2 = Multicomplex{Bool,47,2}(SVector{2}(false, true))
# etc

Multicomplex numbers are often built-up (and mathematically defined) recursively from pairs of lower-order numbers:

m_(n) = a_(n-1) + i_n b_(n-1)

so this route to construction avoids thinking of the number of components (directly, at least). It’s a grey area though, and I agree that there’s definitely times when it’s hard to avoid thinking about the number of components. I do feel it’s generally going to be preferable to expose N to the user rather than C though.

I’ve got some checks in the constructors that I didn’t show to avoid invalid pairs of N and C - and you’d still need similar checks to ensure C was a power of two.

That’s helpful to know about Tuples and arithmetic - it would be good to reduce dependence on other packages. One other area I use StaticArrays is when gradually constructing vectors via an MVector, e.g. to pad a smaller vector with zeros:

v = zeros(MVector{C,T})
for i=1:C0
    v[i] = m.value[i]

I don’t know if there’d be a cleaner way of doing this, particularly with Tuple?

Thanks again - I’m appreciating the discussion!

I can see why that’s more aesthetically pleasing, yes.

As for this:

you could try something like this (I mean, I’m not certain what the interface should be, but this is equivalent, I think):

function padto(t::NTuple{N, T}, ::Val{P}) where {N, P, T}
    m = P - N
    return (t..., ntuple(_->zero(T), m)...)

I was actually not sure if the compiler would be able to optimize away the P - M, but I think it did. So this is if you want to pad to a certain length. If you want to pad with a certain number of zeros, it’s easier:

function padwith(t::NTuple{N, T}, ::Val{P}) where {N, P, T}
    return (t..., ntuple(_->zero(T), P)...)

(You call them as padto(t, Val(7)) to pad to length 7, etc.)

Great - that looks really useful. I’ll definitely have a think about switching across. Thanks for all the help!