Cube root of a negative number?

[shihnh]

> tt <- (-1)^(1/3)
> tt
[1] NaN

Why tt is NaN? For the root in real number, tt should be -1. Did I make any mistake?

[technocrat]

Hi, and welcome!

This post will give you some insight. Basically, -1 is not the same as 0 + 1i

[joels]

-1 has three cube roots (the three solutions to the equation x^3=-1). One of the solutions is -1 and the other two are complex numbers. Going counterclockwise around the complex plane from the positive real axis, the first root is complex, and this is the root that R returns when you work with complex numbers (see below). Since R is not expecting complex numbers in the case of (-1)^(1/3) we need to convert -1 to complex class before taking roots. For example:

as.complex(-1)
#> [1] -1+0i
x = as.complex(-1)^(1/3)   # (-1+0i)^(1/3) will also work
x
#> [1] 0.5+0.8660254i

So we now have one cube root of -1. We can get the other two by rotating x by 120 and 240 degrees (\frac{2}{3} \pi and \frac{4}{3} \pi radians) around the origin in the complex plane.

# Multiplying by the number below will rotate a complex number by 
#  120 degrees (2*pi/3) about the origin 
rot120 = exp(2*pi*(0+1i)/3)  

# Get the other two cube roots of -1 by rotating the first cube root by 
#  120 degrees and 240 degrees about the origin
x = x * rot120^(0:2)

x
#> [1]  0.5+0.8660254i -1.0+0.0000000i  0.5-0.8660254i
x^3
#> [1] -1+0i -1+0i -1+0i

We can also plot the cube roots to see where they're located in the complex plane. R's plot function understands complex numbers and automatically uses the real parts for the x-axis values and the imaginary parts for the y-axis values.

par(pty="s") # square aspect ratio

# Set up plot region
plot(NA, xlim=c(-1,1), ylim=c(-1,1), asp=1, ylab="Im", xlab="Re", las=1)
clip(-1.05,1.05,-1.05,1.05)
abline(h=0, col="grey50")
abline(v=0, col="grey50")

# Plot cubes of cube roots of -1 (all three equal -1)
points(x^3, pch=16, col="red", cex=2)

# Plot cube roots of -1
points(x, pch=21:23, bg="blue")

[shihnh]

Thanks technocrat and joels. The only way i found about cubic root transformation in real domain is to define a function such as

Math.cbrt <- function(x) {
    sign(x) * abs(x)^(1/3)
}
x <- c(-1, -8, -27, -64)
Math.cbrt(x)
# [1] -1 -2 -3 -4

which originally came fromJosh OBrien.
A little bit cumbersome though. It is difficult to generalize to all x^(a), for both x and a are real.

[joels]

What about something like the code below? It returns nth roots of a number. I've set it up to return all roots as well as real roots (which are those roots with Imaginary part equal to 0), just for illustration, but you can modify it to return just the real roots if desired.

# Return the nth roots of a number
# Arguments
# x: A number whose nth roots we want to determine. Can be real or complex.
# n: A positive integer representing the desired nth roots.
real.roots = function(x, n) {
  x = as.complex(x)^(1/n)
  rot.angle = exp(2*pi*(0+1i)/n)
  
  xall = x * rot.angle^(0:(n-1))
  
  list(
    real.roots = Re(xall[purrr::map_lgl(xall, ~isTRUE(all.equal(Im(.x), 0)))]),
    all.roots = xall
  )
}

# Check 
r = real.roots(-1, 3)
r[[2]]^3
#> [1] -1+0i -1+0i -1+0i
r = real.roots(1+2i,4)
r[[2]]^4
#> [1] 1+2i 1+2i 1+2i 1+2i

real.roots(-1, 3)
#> $real.roots
#> [1] -1
#> 
#> $all.roots
#> [1]  0.5+0.8660254i -1.0+0.0000000i  0.5-0.8660254i

real.roots(1+1i, 3)
#> $real.roots
#> numeric(0)
#> 
#> $all.roots
#> [1]  1.0842151+0.2905146i -0.7937005+0.7937005i -0.2905146-1.0842151i

real.roots(1i^1i, 4)
#> $real.roots
#> [1]  0.6752319 -0.6752319
#> 
#> $all.roots
#> [1]  0.6752319+0.0000000i  0.0000000+0.6752319i -0.6752319+0.0000000i
#> [4]  0.0000000-0.6752319i

real.roots(1, 5)
#> $real.roots
#> [1] 1
#> 
#> $all.roots
#> [1]  1.000000+0.0000000i  0.309017+0.9510565i -0.809017+0.5877853i
#> [4] -0.809017-0.5877853i  0.309017-0.9510565i

real.roots(-27, 3)
#> $real.roots
#> [1] -3
#> 
#> $all.roots
#> [1]  1.5+2.598076i -3.0+0.000000i  1.5-2.598076i

real.roots(-1.45 + 0.3i, 3)
#> $real.roots
#> numeric(0)
#> 
#> $all.roots
#> [1]  0.6356515+0.9460756i -1.1371512+0.0774525i  0.5014998-1.0235281i

real.roots(16, 4)
#> $real.roots
#> [1]  2 -2
#> 
#> $all.roots
#> [1]  2+0i  0+2i -2+0i  0-2i

^{Created on 2019-12-30 by the reprex package (v0.3.0)}

[shihnh]

That is Great! joels. Inspired by your code, could we extend the nth roots from positive integer n to any rational number, say n/m? I mean, raised x to the power of m first and then come the code in real.roots? I try it and find that it make no difference between n=3 and n=6/2 etc. Of course, negative rational number still works. thanks a lot!