Custom Complex number type for specific Real type

General Usage
,
1

Hello! I have implemented a Julia interface to a C library for a special number type. The real and complex types are handled internally by the C library, so I have Julia wrapper structs that are basically MyNumber and ComplexMyNumber and have already laboriously implemented all of the operators including with promotion, etc (the C library has special operators when promoting a MyNumber to ComplexMyNumber).

I would like for MyNumber <: Real and ComplexMyNumber <: Number, however I see in base/complex.jl:

struct Complex{T<:Real} <: Number
    re::T
    im::T
end

ComplexMyNumber, which is handled internally by the C library, cannot be defined as two separate MyNumbers for the real and imaginary parts. In my wrapper, ComplexMyNumber and MyNumber are just structs containing Ptrs to the C objects. So, is there some way that I can basically force that the type Base.Complex{MyNumber} = ComplexMyNumber so that whenever a Complex{MyNumber} is defined, it uses my special struct ComplexMyNumber instead?
.

2

Technically you can make a constructor not return the specified type.

Base.Complex{MyNumber}(x::MyNumber) = ComplexMyNumber(x)
Base.Complex{MyNumber}(x::MyNumber, y) = ComplexMyNumber(x, convert(MyNumber, y))

But you shouldn’t have to. And its generally pretty weird and unexpected when people do so

You can overload

Base.complex(x::MyNumber) = ComplexMyNumber(x)
Base.complex(::Type([MyNumber}) = ComplexMyNumber

Which is the function that is supposed to be called to determine what complex type is used for any given real type. (If someone is instead hardcoding it to Compelx{MyNumber} then it is arguably a bug in their code)
From the docs:

complex(T::Type)

Return an appropriate type which can represent a value of type T as a
complex number.

and you have already implemented all
“all of the operators” which is generally what I would suggest doing.
Then anything that accepts Numbers in general should work for your type

What code right now exists that is causing you problems that you want to solve in this way?

2 Likes
3

I suppose it was just an instinct to subtype, alias, or otherwise relate ComplexMyNumber to the existing Complex like you would have MyNumber subtype Real, but there’s really no good reason for it. Complex is a parametric UnionAll type that contains concrete DataType types with specified parameters. Since ComplexMyNumber cannot have the structure of Complex{MyNumber}, they cannot be the same concrete type, they cannot share methods, and Complex cannot contain ComplexMyNumber. In fact, I would suggest making the Complex{MyNumber} constructor throw a descriptive error if you should always use ComplexMyNumber instead. There’s no AbstractComplex type for complex numbers like the abstract DataType type Real for real numbers, which is why Complex subtypes Number; that does also mean the methods annotated by Complex in complex.jl is specific to that struct, and you’ll need to make analogous methods for ComplexMyNumber.

4

Let’s say we have this definition of ComplexMyNumber.

julia> struct ComplexMyNumber <: Number
           ptr::Ptr{Float64}
       end

We can override Base.getproperty for this type.

julia> function Base.getproperty(c::ComplexMyNumber, s::Symbol) 
           if s == :re
               unsafe_load(getfield(c, :ptr), 1)
           elseif s == :im
               unsafe_load(getfield(c, :ptr), 2)
           else
               getfield(c, s)
           end
       end

Then we could do the following to create virtual properties that do not correspond to actual fields.

julia> v = Float64[3,4]
2-element Vector{Float64}:
 3.0
 4.0

julia> cmn = ComplexMyNumber(pointer(v))
ComplexMyNumber(Ptr{Float64} @0x00007f42a82663e0)

julia> cmn.re
3.0

julia> cmn.im
4.0

We can then define conversion from ComplexMyNumber to Complex{Float64}.

julia> Base.convert(::Type{Complex{Float64}}, c::ComplexMyNumber) =
           Complex{Float64}(c.re, c.im)

julia> convert(Complex{Float64}, cmn)
3.0 + 4.0im

We can then take advantage of implicit conversion such as when pushing or assigning into an array.

julia> complex_vector = Complex{Float64}[]
ComplexF64[]

julia> push!(complex_vector, cmn)
1-element Vector{ComplexF64}:
 3.0 + 4.0im

julia> v2 = Float64[9,10]
2-element Vector{Float64}:
  9.0
 10.0

julia> cmn2 = ComplexMyNumber(pointer(v2))
ComplexMyNumber(Ptr{Float64} @0x00007f42a7fde930)

julia> cmn[1]
ComplexMyNumber(Ptr{Float64} @0x00007f42a82663e0)

julia> complex_vector
1-element Vector{ComplexF64}:
 3.0 + 4.0im

julia> complex_vector[1]
3.0 + 4.0im

julia> complex_vector[1] = cmn2
ComplexMyNumber(Ptr{Float64} @0x00007f42a7fde930)

julia> complex_vector[1]
9.0 + 10.0im
1 Like
5

There isn’t actually any “problem” in my code right now, it’s just that MyNumber and ComplexMyNumber do not really conform with the Complex{T<:Real} defined in Base. And so I’m not sure about saying MyNumber <: Real, when MyNumber cannot satisfy everything that a Real type can in Base. However in theory MyNumber and ComplexMyNumber are “Real” and “Number” types respectively, and for usability of the package they should be allowed where these abstract types are allowed, but I guess that does mean I need to play some tricks on the back end to prevent undefined behavior.

This feels like the safest solution, since ComplexMyNumber is indeed not a Complex{MyNumber} as you say, and it prevents users from attempting to define a ComplexMyNumber in this way. But does this cover all bases when defining MyNumber <: Real and ComplexMyNumber <: Number?

Though the implementation of MyNumber is much more complicated than this on the C end, I think I could make something like this work for ComplexMyNumber. However because Complex{MyNumber} is not actually a ComplexMyNumber, I do have some concerns about undefined behavior from operations defined in Base with this method

6

I’m not exactly sure what covers all bases means here, but I would say no because you can see many fundamental complex number methods like imag, real, (these 2 I would put in an nonexistent AbstractComplex interface) abs, angle, conj, etc. that all dispatch on Complex, which won’t apply to ComplexMyNumber. You would basically have to reimplement complex.jl at minimum, then any method that dispatches on Complex instead of Number. If it is even possible, converting between Complex and ComplexMyNumber would add some processing time but you can leverage the existing Complex methods. I’m not sold on implementing getproperty to emulate Complex’s .re and .im fields because the real and imag functions are used instead to be consistent with methods for other Numbers; check complex.jl, the only time .re and .im shows up instead is a show method with Complex{Bool}, and I’m not convinced that wasn’t a typo.

7

I mean, at what point am I allowed to say that something is a Real and something is a Number?

All of these methods you mention (real, imag, angle, abs, conj) I have implemented already ( the C library provides these methods). The library also comes with some methods even not in Base (polar, rect, sinhc). I just looked through complex.jl and most of these methods I have already implemented, and the few ones I haven’t yet I can do so.

8

I would definitely favor finishing wrapping as much of the library as possible in base Julia functions and possibly new specific functions for consistent results, and include the Complex conversion as a pragmatic backup. I’m guessing the conversion can also help performance by storing results, specifiically I’m guessing ComplexMyNumber has its structure because computing the real and imaginary part is not cheap enough to do immediately upon instantiation.

Roughly distinguished by whether it’s representing a real number, more precisely by their methods (for example, you want real(::Real) to be like identity). Look at the tree diagram in Numbers · The Julia Language, all base Julia numerical types are Real except for Complex which is a composite type parameterized by <:Real. Don’t worry about type hierarchies not reflecting mathematical sets perfectly, types are for instantiation and dispatch first and foremost.

1 Like
9

What you really want here is the proposed AbstractComplex in this issue: abstract complex? · Issue #33246 · JuliaLang/julia · GitHub. Unmerged implementation here: https://github.com/JuliaLang/julia/pull/35587.

3 Likes
10

Just to make sure I understand correctly, you recommend allowing Complex{MyNumber} as a backup? This in theory should work for all the Base methods I overloaded, but for the C library specific methods not in Base, would not be allowed. It certainly would be slower to use this backup.

This seems great! Only question I’d have is how to make sure that I can still declare MyNumber <: Real and ensure that the AbstractComplex representation is used by default instead of Complex{MyNumber} (preventing/making it very difficult for users to say Complex{MyNumber} over ComplexMyNumber)

For now at least I can call ComplexMyNumber <: Number and just wrap as much of the library as possible in Julia base functions, and decide how to handle or just prevent Complex{MyNumber}