cbrt(::Complex)

[stevengj]

As noted on discourse, this method is missing.

An obvious (but inconsistent, see below) implementation would be:

cbrt(z::Complex) = cbrt(abs(z)) * cis(angle(z)/3)

but maybe there is a faster way?

Might also be worth having a fallback method too: (this is inconsistent too)

cbrt(x::Number) = x^(1//3)

[stevengj]

Hmm, one issue with the implementation above is that it is inconsistent with the cbrt for real functions, since cbrt(negative real) is a negative real value, but cbrt(::Complex) using the definition above would be complex:

julia> Base.cbrt(z::Complex) = cbrt(abs(z)) * cis(angle(z)/3)

julia> cbrt(-3)
-1.4422495703074083

julia> cbrt(-3+0im)
0.7211247851537043 + 1.2490247664834064im

julia> (-3+0im)^(1//3)
0.7211247851537042 + 1.2490247664834064im

It's not obvious to me what definition we would want for cbrt(::Complex), but it needs to coincide with cbrt(::Real) for real values.

[antoine-levitt]

Having this consistent with real input will necessarily put a branch cut at a non-standard place (different than R- used for the other usual functions). If you compute the cube root of a complex number you will most likely want to specify your branch cut explicitly, so I don't think this function should be defined at all.

[StevenSiew]

In Mathematica, cbrt is only defined to return the real value of cuberoot.

I don't think it is a good idea to define cbrt for complex value because it would cause discontinuity (depending on whether the arguement is Real or Complex). But certainly the documentation for crbt should be updated to mention that

cbrt(-27) == -3 while (-27+0im)^(1/3) == 1.5 + 2.598076211353316im

so cbrt will only return a real result (and never a complex result)

[ettersi]

I'd say the most sensible definition of cbrt(::Complex) is to choose the branch cut [-Inf,Inf] * im \ {0} and fold the left half-plane towards the negative real axis and the right half-plane towards the positive real axis.

Here's an implementation based on this definition.

function Base.cbrt(z::Complex)
    if isreal(z)
        return Complex(cbrt(real(z)),0*imag(z))
    elseif isinf(real(z)) && !isinf(imag(z))
        return Complex(real(z),0*imag(z))
    else
        if !signbit(real(z))
            return cbrt(abs(z)) * cis(angle(z)/3)
        else
            T = float(real(typeof(z)))
            return cbrt(abs(z)) * cis(angle(z)/3 + copysign(2*T(π)/3,imag(z)))
        end
    end
end

Tests, imitating the tests for sqrt(::Complex) from test/complex.jl, lines 211 - 250:

@testset "cbrt" begin
    @testset for T in (Float16,Float32,Float64,BigFloat,Int8,Int16,Int32,Int64,Int128,BigInt)
        @inferred cbrt(zero(Complex{T}))
    end

    @testset for T in (Float32, Float64, BigFloat)
        x = Complex{T}(1//3 + 1//4*im)
        @test cbrt(x)^3 ≈ x
    end

    for x = ( ComplexF64(s1//3 + s2//4*im) for s1 in (-1,1), s2 in (-1,1) )
        @test all(signbit.(reim(cbrt(x))) .== signbit.(reim(x)))
    end

    @test isequal(cbrt(complex( 0.0, 0.0)), complex( 0.0, 0.0))
    @test isequal(cbrt(complex( 0.0,-0.0)), complex( 0.0,-0.0))
    @test isequal(cbrt(complex( 0.0, Inf)), complex( Inf, Inf))
    @test isequal(cbrt(complex( 0.0,-Inf)), complex( Inf,-Inf))
    @test isequal(cbrt(complex( 0.0, NaN)), complex( NaN, NaN))

    @test isequal(cbrt(complex(-0.0, 0.0)), complex(-0.0, 0.0))
    @test isequal(cbrt(complex(-0.0,-0.0)), complex(-0.0,-0.0))
    @test isequal(cbrt(complex(-0.0, Inf)), complex(-Inf, Inf))
    @test isequal(cbrt(complex(-0.0,-Inf)), complex(-Inf,-Inf))
    @test isequal(cbrt(complex(-0.0, NaN)), complex( NaN, NaN))

    @test isequal(cbrt(complex( 5.0, 0.0)), complex( cbrt(5.0), 0.0))
    @test isequal(cbrt(complex( 5.0,-0.0)), complex( cbrt(5.0),-0.0))
    @test isequal(cbrt(complex(-5.0, 0.0)), complex(-cbrt(5.0), 0.0))
    @test isequal(cbrt(complex(-5.0,-0.0)), complex(-cbrt(5.0),-0.0))

    @test isequal(cbrt(complex( Inf, 0.0)), complex( Inf, 0.0))
    @test isequal(cbrt(complex( Inf,-0.0)), complex( Inf,-0.0))
    @test isequal(cbrt(complex( Inf, 5.0)), complex( Inf, 0.0))
    @test isequal(cbrt(complex( Inf,-5.0)), complex( Inf,-0.0))
    @test isequal(cbrt(complex( Inf, Inf)), complex( Inf, Inf))
    @test isequal(cbrt(complex( Inf,-Inf)), complex( Inf,-Inf))
    @test isequal(cbrt(complex( Inf, NaN)), complex( Inf, NaN))

    @test isequal(cbrt(complex(-Inf, 0.0)), complex(-Inf, 0.0))
    @test isequal(cbrt(complex(-Inf,-0.0)), complex(-Inf,-0.0))
    @test isequal(cbrt(complex(-Inf, 5.0)), complex(-Inf, 0.0))
    @test isequal(cbrt(complex(-Inf,-5.0)), complex(-Inf,-0.0))
    @test isequal(cbrt(complex(-Inf, Inf)), complex(-Inf, Inf))
    @test isequal(cbrt(complex(-Inf,-Inf)), complex(-Inf,-Inf))
    @test isequal(cbrt(complex(-Inf, NaN)), complex(-Inf, NaN))

    @test isequal(cbrt(complex( NaN, 0.0)), complex( NaN, 0.0))
    @test isequal(cbrt(complex( NaN,-0.0)), complex( NaN,-0.0))
    @test isequal(cbrt(complex( NaN, Inf)), complex( NaN, NaN))
    @test isequal(cbrt(complex( NaN,-Inf)), complex( NaN,-NaN))
end

[ettersi]

Any chance this will get approved if I turn the above into a PR?

[stevengj]

I think that we should be cautious about implementing this, since the "correct" choice is non-obvious and the difference from z^(1/3) could be extremely confusing, until we have specific applications to shed light on what is actually needed (if anything) here.

[stevengj]

I've been thinking about this some more in the context of #47565, and it seems like there is some value in defining a cbrt(z::Complex) function that goes to cbrt(x::Real) as z ⟶ x. In particular, this is (I think?) equivalent to @ettersi's definition above:

cbrt(z::Complex) = signbit(real(z)) ? -(-z)^(1//3) : z^(1//3)

For example, this would allow one to define a real cube root of a diagonalizable real matrix A=XΛX⁻¹ (#47513) as B = X cbrt.(Λ) X⁻¹. Importantly, this definition is consistent with what you would get for Hermitian matrices by taking cbrt of the real eigenvalues. (You wouldn't want to compute it this way because X might be badly conditioned, but at least this gives a definition, which even extends in principle to non-diagonalizable matrices by taking the the limit.)

[theWiseAman]

Hmm, one issue with the implementation above is that it is inconsistent with the cbrt for real functions, since cbrt(negative real) is a negative real value, but cbrt(::Complex) using the definition above would be complex:

julia> Base.cbrt(z::Complex) = cbrt(abs(z)) * cis(angle(z)/3)

julia> cbrt(-3)
-1.4422495703074083

julia> cbrt(-3+0im)
0.7211247851537043 + 1.2490247664834064im

julia> (-3+0im)^(1//3)
0.7211247851537042 + 1.2490247664834064im

It's not obvious to me what definition we would want for cbrt(::Complex), but it needs to coincide with cbrt(::Real) for real values.

Actually there would be three answers (3 cube roots in the complex domain) in this case one of which would lie along the real axis.

[theWiseAman]

I've been thinking about this some more in the context of #47565, and it seems like there is some value in defining a cbrt(z::Complex) function that goes to cbrt(x::Real) as z ⟶ x. In particular, this is (I think?) equivalent to @ettersi's definition above:

cbrt(z::Complex) = signbit(real(z)) ? -(-z)^(1//3) : z^(1//3)

For example, this would allow one to define a real cube root of a diagonalizable real matrix A=XΛX⁻¹ (#47513) as B = X cbrt.(Λ) X⁻¹. Importantly, this definition is consistent with what you would get for Hermitian matrices by taking cbrt of the real eigenvalues. (You wouldn't want to compute it this way because X might be badly conditioned, but at least this gives a definition, which even extends in principle to non-diagonalizable matrices by taking the the limit.)

@stevengj Please check my Draft PR (#47860). I think to make nthroot(z::Complex) consistent when x::Real we can throw real value from the array of roots. There is sure to be at least one real value when x < 0 && n::odd else if x < 0 && n::even we need to throw error that nthRoot falls in Complex Domain as implemented in sqrt(). Also, I don't think we need cbrt() if nthroot() is completed.

[stevengj]

@theWiseAman, we don't want an array of roots. The whole point of cbrt is to choose a particular root which is different from the usual principal root for negative arguments.

[MilesCranmer]

Ping on this. Ran into this today

[stevengj]

@MilesCranmer, what behavior did you want cbrt to have for complex numbers in your application?

[MilesCranmer]

I won’t need a particular branch cut for my work as it would be for SymbolicRegression.jl where the search could learn an additional rotation in the complex plane if needed. (i.e., any branch cut would do)

I suppose my vote would be for matching the same branch cut as x ^ (1/3) as that’s what I had assumed initially would be true.

[stevengj]

I suppose my vote would be for matching the same branch cut as x ^ (1/3) as that’s what I had assumed initially would be true.

That choice seems out of the question for me, since it would be inconsistent with cbrt(::Real). i.e. you would get totally different resultx for cbrt(x) and cbrt(complex(x)).

But then, having cbrt(z::Complex) behave differently from z^(1/3) is also surprising. So, that argues in favor of continuing to not implement this method.

[MilesCranmer]

Yeah it actually looks like this – a lack of any implementation – seems to be the standard:

• Python doesn't implement it: https://bugs.python.org/issue44357 (related discussion)
  • I also checked scipy.special which seems to skip it.
• C++ stdlib implements a complex sqrt: https://en.cppreference.com/w/c/numeric/complex/csqrt but does not implement a complex cbrt: