Julia: 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.