Add abstract complex

[MasonProtter]

Closes #33246

I said I'd do this a while ago and next got around to it. I'm sure there's a lot more machinery that can be put into this, but I figured this was approximately the minimal set of methods that AbstractComplex should satisfy, with the rest being up to whoever implements other complex types.

[MasonProtter]

Something perhaps worth discussing: should we have Real <: AbstractComplex? I think the answer is no, but I do think that there could be some advantage to having an abstract type for commutative numbers whereas the rest of Number's subtypes are assumed to be non-commutative.

This slight possible advantage is likely not worth it though. Interested to hear if anyone has thoughts or feelings on the matter.

[kimikage]

From a programming perspective rather than a mathematical perspective, why is AbstractComplex not parametric?

{T} is not good semantically, but I think it is a necessary evil in order to make AbstractComplex more widely used. Do I have to use, for example, real(::Type{<:AbstractComplex}) to support AbstractComplex? What a wonderful abstraction it is!

[jonas-schulze]

As real numbers are a strict subset of complex numbers, I’d be in favor of Real <: AbstractComplex.

[StefanKarpinski]

Regarding the type parameter: yes, it seems like AbstractComplex should have a type parameter: it tells you what kind of value you get when you ask for things like the real or imaginary component of a complex value or ask for its absolute value.

We're not doing Real <: AbstractComplex. I've explained this before in various places: while the real numbers are identified with a subset of the complex numbers, they are not actually a subset of the complex numbers—they are (field) isomorphic to a subset. In the standard definition, the complex numbers are constructed as a pair of reals, so you need to have the reals before you can even define the complex numbers, so logically the reals cannot be a subset of the complexes. Yes, there are other non-constructive ways to define the complex numbers that skip right over the reals, but on a computer constructive approaches are the only option. Trying to make something a subset of the a thing that it is defined in terms of does not lead anywhere good. The right way to identify the reals with a subset of complex numbers computationally is to make conversion work correctly and make them behave the same.

[MasonProtter]

Regarding the type parameter: yes, it seems like AbstractComplex should have a type parameter: it tells you what kind of value you get when you ask for things like the real or imaginary component of a complex value or ask for its absolute value.

The thing that worried me about the type parameter in AbstractComplex was the case where someone might want to do something like

struct HeterogeneousComplex{T, U}
   re::T
   im::U
end

but I guess this is a rather small corner case and whoever implemented such a thing can just override the Base methods still. I'll add in the type parameter.

[MasonProtter]

@Roger-luo are you still interested in collaborating on this? If you or anyone else is interested, I could use some help coming up with test cases and docs.

I'm also not sure what to do with promotion rules so (if anything) suggestions on that front would be appreciated.

[kimikage]

abstract type AbstractComplex{Tre<:Real, Tim<:Real} <: Number end

struct Complex{T} <: AbstractComplex{T, T}
    re::T
    im::T
end

struct HeterogeneousComplex{T, U} <: AbstractComplex{T, U}
    re::T
    im::U
end

😄

[MasonProtter]

Yes, that is a possibility, but would such a thing be worth the complication it adds? I suspect no but I'm open to arguments.

[kimikage]

I cannot imagine other practical complex types except the polar representation. So, I think one parameter is sufficient. (Since trailing parameters can be omitted, I also think the 2-parameter type is not so "complex".)

[kimikage]

I don't really understand what properties the AbstractComplex should have, i.e. the degree of abstraction. Which document should I read to understand the difference from two years ago? (cf. #26666 (comment))

As an end user, I welcome this change, but as a package developer, I'm reluctant to support AbstractComplex. In practice, what advantage does AbstractComplex have over Union{Complex, PolarComplex}?

[MasonProtter]

In practice, what advantage does AbstractComplex have over Union{Complex, PolarComplex}?

The difference is that if you are using AbstractComplex, your package doesn't need to know about PolarComplex or whatever in order to support it. It can come from some random package you've never heard of, but inter-operate with your code if it satisfies the right interface.

[kimikage]

Yes. So I'm trying to understand "the right interface". If it satisfies the right interface, my code will work with the Number instead of AbstractComplex.:stuck_out_tongue:

I think Complex should be the reference implementation of AbstractComplex. The parameter {T} is based on this idea. So IMO, dual numbers, split-complex numbers and so on should not be included in AbstractComplex, and that should be clearly documented. Otherwise, we will add SuperAbstractComplex in the near future.

I'm not against types like SuperAbstractComplex (apart from the name), but the type hierarchy shouldn't change often.

[MasonProtter]

I think Complex should be the reference implementation of AbstractComplex. The parameter {T} is based on this idea. So IMO, dual numbers, split-complex numbers and so on should not be included in AbstractComplex, and that should be clearly documented. Otherwise, we will add SuperAbstractComplex in the near future.

ForwardDiff.Dual gets by just fine as a subtype of Real. Why shouldn't a WirtingerDual subtype AbstractComplex?

Split complex numbers are a bit of a harrier situation. I think you're be right to be suspicious of them being <:AbstractComplex.

[Roger-luo]

Hi, sorry for this late reply.

My originally motivation for AbstractComplex interface is mainly for complex types that cannot be represented by two real numbers. This is the original idea for this.

A more concrete use case is the symbolic complex type, that treats a symbol as a complex number instead of two variable that is subtype of Real. Thus as far as I can understand it doesn't make sense to have type parameters for AbstractComplex in this case since there is nothing really need to be attached to its parent type. Assuming we have a symbolic complex type

struct SComplex <: AbstractComplex
   ex
end

x = SComplex(:z) 

the return value of real should be real(:z) which is of SReal type that contains an expression, but in symbolic context a symbolic complex number can sometimes also be a constant of other numeric types. This makes the return value not usually very type stable and it's not necessary to be type stable either.

Expr
  head: Symbol call
  args: Array{Any}((2,))
    1: Symbol real
    2: SComplex

thus you can't really know what makes sense to put into AbstractComplex's type parameter.
The other case I could think of is a pure imaginary type Im{T}. I feel AbstractComplex should be something similar to Number and Real, which is just a general separation of the interface instead of being too verbose. And this is more consistent with other abstract number types I believe.

[Roger-luo]

What I'm proposing here is

struct Complex{T} <: AbstractComplex
   re::T
   im::T
end

and for other types

struct SComplex <: AbstractComplex
   ex
end

struct Im{T} <: AbstractComplex
   im::T
end

struct Polar{T} <: AbstractComplex
   re::T
   im::T
end

I don't think we really need a type parameter, even if we need one, the better way is to use typeof(real(x)) directly which covers all dynamic cases in symbolic computation.

[kimikage]

If Complex{T} didn't exist yet, I would have really agreed to the non-parametric AbstractComplex. I am not the beneficiary, so I will wait and see without interfering.:speak_no_evil:

[kalmarek]

yet another (semi-symbolic) representation of complex numbers: In https://github.com/kalmarek/Cyclotomics.jl I use

julia> x = E(5)^2 + E(5)^3
 ζ₅² +ζ₅³

julia> dump(x)
Cyclotomics.Cyclotomic{Int64,SparseArrays.SparseVector{Int64,Int64}}
  n: Int64 5
  coeffs: SparseArrays.SparseVector{Int64,Int64}
    n: Int64 5
    nzind: Array{Int64}((2,)) [3, 4]
    nzval: Array{Int64}((2,)) [1, 1]

i.e. represent sums of roots of unity as elements of the appropriate vector space over T

julia> E(45)^9
 ζ₅

julia> E(45)^5
 -ζ₉⁴ -ζ₉⁷

julia> E(45)^5 + 2E(45)^9
 ζ₄₅² +ζ₄₅⁸ +ζ₄₅¹¹ +ζ₄₅¹⁷ -2*ζ₄₅²⁴ +ζ₄₅²⁶ +ζ₄₅²⁹ +ζ₄₅³⁸ -2*ζ₄₅³⁹ +ζ₄₅⁴⁴

clearly this is not (re, im) representation and my type is parametrized by the type of coefficients + the type of storage:

julia> typeof(1.0*E(5,3))
Cyclotomics.Cyclotomic{Float64,SparseArrays.SparseVector{Float64,Int64}}

I'd like to make it into AbstractComplex{T} though, but this T is very much different from T of Complex{T}...

[Keno]

Bump :) I ran into an application where having a well supported Polar representation type would be useful :). What remains to be decided here?

[KristofferC]

I think showing any usefulness in practice at having the abstract type, cf #33246 (comment)

[Roger-luo]

Bump :) I ran into an application where having a well supported Polar representation type would be useful :). What remains to be decided here?

@Keno I think some people (like me) don't think there should be a type parameter for the abstract type since it can be just a constant or some struct with an empty field like Number, but the implementation here does.

The rest I think looks good to me.

[MasonProtter]

I don't have any strong feelings one way or another on whether or not it should be parameterized. I can see benefits and downsides to either. Keno, since you want this, do you have any intuitions regarding parameterizing AbstractComplex or not?

[StefanKarpinski]

I think some people (like me) don't think there should be a type parameter for the abstract type since it can be just a constant or some struct with an empty field like Number, but the implementation here does.

I think there should be a type parameter: the type parameter tells you what kind of value to expect when you ask for the real or imaginary parts of the number. That does not force the implementation to store a value of that type or any value for that matter. If the value is a constant you still have to return a value when asking for the real/imaginary parts and whatever that type is is what the type parameter should be.

[JeffBezanson]

Right, you can always set the type parameter to Any if you really need to.

I don't understand the SComplex use case. If real(SComplex(z)) gives the symbolic expression real(z), why not skip the SComplex type and just use the expression real(z) to begin with? Or, if you have symbolic expressions, use x + y*im or r*exp(theta*im).

[MasonProtter]

I think the idea is that in your symbolic system, you'd write something like

z = SComplex(:z)
x = SReal(:x)

ex1= simplify(z^2 >= 0)
ex2 = simplify(x^2 >= 0

and have ex1 not be changed but ex2 to just simplify to true. Roger isn't wanting to talk about a complex number as being an object composed of real numbers, but instead as it's own object.

It's directly anaologous to wanting a symbolic representation of an array, rather than an array of symbolic scalars.

[JeffBezanson]

From triage:

• Might as well do this
• We like the type parameter

[Roger-luo]

After some thinking about the parameteric solution, I think actually since there is always a symbolic term type in CAS system, e.g SymbolicUtils has Term, the implementation can be just

struct SComplex <: AbstractComplex{Term}
   ex::Term
end

and Term will take the role to wrap Expr or Number or other things.

[MasonProtter]

So what should the minimal interface be? Just real and imag? I suppose that should be documented.

There's some code here that feels like it could be made more generic, like

julia/base/complex.jl

Lines 290 to 292 in f1915bc

	muladd(z::Complex, w::Complex, x::Complex) =
	Complex(muladd(real(z), real(w), -muladd(imag(z), imag(w), -real(x))),
	muladd(real(z), imag(w), muladd(imag(z), real(w), imag(x))))

But I'm not sure what the generalization would be. Would the fallbacks just convert to Complex? I feel like that'd be annoying for some people implementing their own complex numbers, but perhaps convenient for others.

Perhaps we could make part of the interface that your complex type needs a way to take in real and imaginary parts to make an instance, e.g. every MyComplex{T} <: AbstractComplex should have a constructor

MyComplex{T}(; real, imag)

Then for instance, muladd could be written

function muladd(z::AbstractComplex, w::AbstractComplex, x::AbstractComplex)
    T = promote_type(typeof(z), typeof(w), typeof(x))
    T(real = muladd(real(z), real(w), -muladd(imag(z), imag(w), -real(x))),
      imag = muladd(real(z), imag(w),  muladd(imag(z), real(w),  imag(x))))
end

Users could then make their own more efficient methods if they wanted.

[Roger-luo]

Perhaps we could make part of the interface that your complex type needs a way to take in real and imaginary parts to make an instance, e.g. every MyComplex{T} <: AbstractComplex should have a constructor

I think so, this is consistent with other things in Base too.

[MasonProtter]

Requiring people to make such a constructor with kwargs seemed like a nonstandard thing to do, so instead I'm going to try something like

for unary in (:conj, :inv, :+, :-, :sqrt, :cis, :log, :cispi, :angle, :log10, :log2, :exp, :expm1, :log1p, :exp2, :exp10,
              :sin, :cos, :tan, :asin, :acos, :atan, :sinh, :cosh, :tanh, :asinh, :acosh, :atanh, :float, :big)
    @eval $unary(z::AbstractComplex) = (typeof(z) ∘ $unary ∘ Complex)(z)
end

for binary in (:+, :-, :*, :/, :^, :muladd)
    @eval $binary(z::AbstractComplex, w::AbstractComplex) = (promote_typeof(z, w) ∘ $binary)(Complex(z), Complex(w))
    @eval $binary(z::Real, w::AbstractComplex) = (promote_typeof(z, w) ∘ $binary)(z, Complex(w))
    @eval $binary(z::AbstractComplex, w::Real) = (promote_typeof(z, w) ∘ $binary)(Complex(z), w)
end

This way, to get a generic implementation, you just need to be able to convert back and forth between your type and Complex.

[oscardssmith]

bump. What is needed to get this merged?

[daanhb]

While we're at it, is it an idea to implement an iscomplex function, analogously to isreal?

[oscardssmith]

what behavior would it have?

[daanhb]

It's only tangential to this PR, sorry, but I'll explain where I'm coming from. I was trying to define the sets of natural numbers, of integers, real numbers and rational numbers in DomainSets.jl (JuliaApproximation/DomainSets.jl#95). Base provides isinteger and isreal functions, which return true if the argument numerically equals an integer or a real number respectively, regardless of its type. For example, isinteger(5.0) is true, so I can do convert(Int, 5.0) with success.

But there are no isrational and iscomplex functions. These would also be less useful, because there is only one Rational type and one Complex type. That might change if there are multiple complex types, hence my suggestion to add them here.

For a number, just like isreal, I'd expect iscomplex to be true if the argument numerically equals a point in the complex plane. Again like isreal, but unlike isinteger, it should probably also be true for arrays of such numbers. Also, iscomplex(x) being true signals to me that I can take the real and imaginary parts of x, and in that case I would expect isreal(real(x)) and isreal(imag(x)) to be true.

Still might not be used much :-)

[JeffreySarnoff]

what behavior would it [iscomplex] have?
it could be

iscomplex(x::Number) = false
iscomplex(x::AbstractComplex) = true