Cube root of a real number

[ivan-pi]

Related to #150 (non-special mathematical functions)

cbrt - Cube root of a real number

Description

Returns the cube root of the real number (x), that is a number (y) such that (y^3 = x).

Syntax

y = cbrt(x)

Arguments

x: A real number (x).

Return value

Returns the value (\sqrt[3]{x}), the result is of the type real and has the same kind as x.

Example

program demo_cbrt
use stdlib_experimental_math, only: cbrt
print *, cbrt(27.), cbrt(-27.)  ! outputs 3, -3
end program

---

As seen from the discussion on Discourse this function is semantically more accurate than writing x**(1./3.) which only works for positive real numbers and returns NaN otherwise.

A possible extension would be to allow complex arguments, the return value would then be an array with 3 elements for the 3 cube roots (if the number is real and non-zero, there is one real root and a conjugate pair of complex roots; a complex non-zero value will have three distinct cube roots)

[ivan-pi]

I've already got an implementation ready based upon the NSWC Mathematical Library version and using Fypp for templating of different real kinds.

[jvdp1]

Thanks.

A possible extension would be to allow complex arguments, the return value would then be an array with 3 elements for the 3 cube roots (if the number is real and non-zero, there is one real root and a conjugate pair of complex roots; a complex non-zero value will have three distinct cube roots)

I don't use usually complex numbers, neither cube roots. However, to be in agreement with the intrinsic function sqrt, I think it should accept complex arguments. In this case, is a function still possible, or should it be a subroutine?

[ivan-pi]

Personally, I also don't have a foreseeable usage for complex cube roots.

It should however, be possible to have the the same interface for real and complex cube roots:

interface cbrt
function cbrt_sp(x) result(cbrt)
  complex(sp), intent(in) :: x
  complex(sp) :: cbrt
end function
function cbrt_complex_sp(x) result(cbrt)
  complex(sp), intent(in) :: x
  complex(sp) :: cbrt(3)
end function
end interface

A better option might be to follow the IMSL library CBRT(X) function:

For complex arguments, the branch cut for the cube root is taken along the negative real axis. The argument of the result, therefore, is greater than –π/3 and less than or equal to π/3. The other two roots are obtained by rotating the principal root by 3 π/3 and π/3.

Lahey/Fujitsu Fortran also provides a CBRT function, which returns a single number for real or complex variables.

In a thread on the Intel Fortran Forum, @sblionel shows the compiler does in fact call a special cube root routine for a**(1./3.), but this only works for a positive argument. This makes me wonder whether a simple implementation for real numbers could be simply:

function cbrt_v1(x) result(cbrt)
  real, intent(in) :: x
  real :: cbrt
  if (x >= 0.0) then
    cbrt = x**(1.0/3.0)
  else
    cbrt = -((-x)**(1.0/3.0))
  end if
end function

Digging back further, the behavior of the a**(1./3.) has changed on occasion in the Intel Fortran compiler:

• Change in floating point rounding between Versions 11 and 12 of Fortran compiler
• Why does Version 19 change how cube root is calculated?

Also the gfortran mailing list contains an interesting discussion on cbrt:

• Gfortran and using C99 cbrt for X ** (1./3.)

@certik Do you have an opinion on this one? My understanding is that the goal of providing a CBRT function is to be "mathematically" correct, in the sense that it can also accept negative arguments. The problem however is can the behavior of CBRT differ from X**(1./3.) in terms of representation/accuracy/speed.

[certik]

See my comment on precisely this issue (and the subsequent discussion there): symengine/symengine#1644 (comment)

We should do whatever is the most consistent and document it well. It could be that we need two cbrt functions for the two most common conventions.

[nshaffer]

How about an optional argument to decide which branch you want?

f = cbrt(z) always has 0 <= arg(f) < 2*pi/3
f = cbrt(z, k) has 2*pi*k/3 <= arg(f) < 2*pi*(k+1)/3

[nshaffer]

As for negative real arguments, I am inclined to have cbrt(-x) == -cbrt(x) (for positive x).

Upshot: this is what we learn in school before getting introduced to complex numbers. It's what most people would expect when working with reals only.

Drawback: introduces some "gotcha" potential for people who use mix reals and complex numbers without reading the docs.

[ivan-pi]

How about an optional argument to decide which branch you want?

f = cbrt(z) always has 0 <= arg(f) < 2*pi/3
f = cbrt(z, k) has 2*pi*k/3 <= arg(f) < 2*pi*(k+1)/3

With an elemental procedure, you could pass an array of k's if you wanted different branches. Not sure, if it is useful in practice.

[certik]

The c++ cbrt has the property of cbrt(-x) == -cbrt(x) for positive x, which makes it not compatible with complex functions. It is equivalent to the Surd function in Mathematica, and specifically to the CubeRoot function. So one option is to simply have cbrt to do the same thing, returning a real number.

And if you want the complex number, you can just use **, as in cmplx(x, dp)**(1._dp/3) where x is some real (or complex) variable.

[nshaffer]

Or better yet, cast it to complex first. I have in mind:

cbrt(-8.) == -2.
cbrt(cmplx(-8.)) == (1., sqrt(3.))
cbrt(cmplx(-8.), k=1) == (-2., 0.)

[jvdp1]

The c++ cbrt has the property of cbrt(-x) == -cbrt(x) for positive x, which makes it not compatible with complex functions. It is equivalent to the Surd function in Mathematica, and specifically to the CubeRoot function. So one option is to simply have cbrt to do the same thing, returning a real number.

And if you want the complex number, you can just use **, as in cmplx(x, dp)**(1._dp/3) where x is some real (or complex) variable.

This is also the behavior of Octave cbrt

[ivan-pi]

MATLAB on the other hand does not provide a cbrt function (see How do you enter the command for a cube root?). It does however have the following behavior:

>> (8)^(1/3)

ans =

     2

>> (-8)^(1/3)

ans =

   1.0000 + 1.7321i

>> nthroot(-8,3)

ans =

    -2

>> help nthroot
 nthroot Real n-th root of real numbers.
 
    nthroot(X, N) returns the real Nth root of the elements of X.
    Both X and N must be real, and if X is negative, N must be an odd integer.
 
    Class support for inputs X, N:
       float: double, single

[certik]

It seems Matlab's ^ operator behaves like Fortran's ** operator for complex numbers. The nthroot seems to be like like Mathematica's Surd.

[ivan-pi]

Let's keep the ball rolling. The way I see it now, there are essentially two choices:

1. Simple version (no branch argument) - Code example
2. With optional branch variable - Code example

Comments:

• For a complex root of a real variable the cbrt is more compact than using the power operator

y = cbrt(cmplx(x))
y = cmplx(x)**(1._sp/3)

y = cbrt(z)
y = z**(1._sp/3)

• If we go with the 1st option the user can always retrieve the other two complex roots by rotating the principal root by 3π/3 and π/3. This can be documented with an example:

r = cmplx(-1.,sqrt(3.))/2.
j = cmplx(0.,1.)

z = -8 + 0*j

y1 = cbrt(z)       !  1.0000 + 1.7321i 
y2 = y1 * r        ! -2.0000 - 0.0000i
y3 = y1 * conjg(r) !  1.0000 - 1.7321i 

• With the second option, users who know exactly what they want, can easily recover a complex root from a either real or complex number, e.g.

  write(*,'(A)')  "real to real"
  
  write(*,fmtr) "cbrt( 8._sp) = ", cbrt(8._sp)
  write(*,fmtr) "cbrt(-8._sp) = ", cbrt(-8._sp)

  write(*,'(/,A)')  "real to complex"

  write(*,fmtc) "cbrt(8._sp,k=0) = ", cbrt(8._sp,k=0)
  write(*,fmtc) "cbrt(8._sp,k=1) = ", cbrt(8._sp,k=1)
  write(*,fmtc) "cbrt(8._sp,k=2) = ", cbrt(8._sp,k=2)

  write(*,'(/,A)')  "complex to complex"

  z = cmplx(-8._sp)
  write(*,fmtc) "z = ", z
  write(*,fmtc) "cbrt(z) = ", cbrt(z)
  write(*,fmtc) "cbrt(z,k=0) = ",cbrt(z,k=0)
  write(*,fmtc) "cbrt(z,k=1) = ",cbrt(z,k=1)
  write(*,fmtc) "cbrt(z,k=2) = ",cbrt(z,k=2)
  write(*,fmtc) "cbrt(z,k=3) = ",cbrt(z,k=3)

produces the following output:

real to real
cbrt( 8._sp) =  2.0000
cbrt(-8._sp) = -2.0000

real to complex
cbrt(8._sp,k=0) =  2.0000+0.0000j
cbrt(8._sp,k=1) = -1.0000+1.7321j
cbrt(8._sp,k=2) = -1.0000-1.7321j

complex to complex
z = -8.0000+0.0000j
cbrt(z) =  1.0000+1.7321j
cbrt(z,k=0) =  1.0000+1.7321j
cbrt(z,k=1) = -2.0000-0.0000j
cbrt(z,k=2) =  1.0000-1.7321j
cbrt(z,k=3) =  1.0000+1.7321j