A negative base raised to a rational exponent should not produce a co…

[Geropy]

…mplex number if the denominator of the exponent is odd; the result should be real.

For example -8 ^ (1/3) should return -2, not a complex number

Added is_even() method to Integer to help accomplish this

[isuruf]

I agree that there is a discrepancy between cbrt(-8)==-2 and cbrt(-8.0) == 1.0+sqrt(3.0)*I.
You are suggesting cbrt(-8.0)==2.0, but I'm not sure we want to do that. For example in SymPy,
cbrt(-8) == 2 * (-1)**(1/3) and cbrt(-8).n() == 1.0+sqrt(3.0)*I

@certik, thoughts?

[Geropy]

The result of std::cbrt(-8.0) is -2.0, so I would certainly argue that this would be the case here as well.
I would agree that pow(-8.0, 1.0/3.0) == 1.0+sqrt(3.0)*I just as in the std, but the cbrt function shouldn't be limited by floating-point precision in the exponent

[Geropy]

The result of std::cbrt(-8.0) is -2.0, so I would certainly argue that this would be the case here as well.
I would agree that pow(-8.0, 1.0/3.0) == 1.0+sqrt(3.0)*I just as in the std, but the cbrt function shouldn't be limited by floating-point precision in the exponent

[Geropy]

closed by accident

[isuruf]

ping @certik

[certik]

The cbrt function is currently defined as:

RCP<const Basic> cbrt(RCP<const Basic> &arg)
{
    return pow(arg, div(one, i3));
}

Which makes sense.

Now let's do the calculation:

cbrt(-8) = (-8)**(1/3) = exp(1/3 * log(-8)) = exp(1/3 * (log(8) + I*pi)) = exp(log(8)/3) * exp(I*pi/3) =
         = 2 * (1/2 + sqrt(3)/2 * I) = 1 + sqrt(3)*I

So mathematically, it must be cbrt(-8) = 1 + sqrt(3)*I if we follow the usual convention for complex branch cuts, see e.g.:

https://www.theoretical-physics.com/dev/math/complex.html

I understand that (-2)^3 = -8. So is (1 + sqrt(3)*I)^3 = -8. However, the cube root, if defined as x^(1/3), must return 1 + sqrt(3)*I and not -2.

@Geropy, @isuruf what is the exact definition of cbrt such that cbrt(-8) = -2?

Update: I see how: if we define cbrt(x) = x^(1/3), but we define all branch cuts differently, namely:

log(-8) = log(8) + I*pi + 2*pi*I*k

where k=0 gives the cbrt(-8) = 1 + sqrt(3)*I result, and k=1 gives the cbrt(-8) = -2 result:

cbrt(-8) = (-8)**(1/3) = exp(1/3 * log(-8)) = exp(1/3 * (log(8) + I*pi + 2*pi*I*k)) = exp(log(8)/3) * exp(I*pi*(1+2*k)/3) =
         = 2 * exp(I*pi*(1+2*k)/3)

The last exponential:

exp(I*pi*(1+2*k)/3) 
        = 1/2 + sqrt(3)/2 * I     ... for k=0
        = -1                      ... for k=1
        = 1/2 - sqrt(3)/2 * I     ... for k=2

So k=1 gives the result you want. However, that means that every other complex function in SymEngine will become inconsistent. One would have to essentially redo all derivations in:

https://www.theoretical-physics.com/dev/math/complex.html

But instead of defining:

arg(z) = atan2(Im z, Re z)

one would have to define:

arg(z) = atan2(Im z, Re z) + 2*pi

And one could do this, and then indeed (-8)^(1/3) = -2. But I don't think we should do that. Rather I suggest we do exactly what C, C++, Fortran, Python and any other languages does, which is the
arg(z) = atan2(Im z, Re z) definition, and then one must get (-8)^(1/3) = 1 + sqrt(3)*I.

[certik]

So for consistency, I would keep the current (complex) definitions of sqrt and cbrt in SymEngine.

However, looking at the std::cbrt documentation, it looks like they define cbrt for negative values simply as follows: cbrt(x) = - cbrt(-x) for negative x. That is not consistent with complex numbers (even in C++), but one can do it just for this particular function.

@Geropy what is your application of this? We could provide a new function cbrt_real that only works for real numbers, and works just like in C++. Essentially the std::cbrt is using a different branch cut than the rest of C++ complex functions, in order to make std::cbrt always real for both positive and negative values.

[Geropy]

The cbrt_real function would work fine for me; something that only works for real numbers, but returns a real solution.

Looking into it a bit more, if you want to be more general, you could potentially implement something like the Surd function from Wolfram

https://mathworld.wolfram.com/nthRoot.html
https://reference.wolfram.com/language/ref/Surd.html

[certik]

@Geropy I see. I think you found it. What you want is the Surd function. In the definition, they are very clear, that for positive real x, Surd is identical to x^(1/n). But for other negative values it gives you the real result. Regarding cbrt, Mathematica has https://reference.wolfram.com/language/ref/CubeRoot.html, that gives the real-valued result.

What does "real valued" sqrt(-1) return using Surd?

[Geropy]

In the "Possible issues" of the Surd definition, it says that it is undefined for even n's and negative real x

[angelog0]

I was working on this argument and I came across this discussion.

This is what I did. It should be self explanatory by means of comments. Maybe you find bugs there or how to do things better.

The FIKE algorithm should do better results.

!
! gfortran -std=f2018 -O2 -Wall cbrt_test.f90 -o cbrt_test.out
!
!
! For a discussion of CBRT function implementation see
!
!   https://github.com/symengine/symengine/pull/1644
!
! and this comment
!
!   https://github.com/symengine/symengine/pull/1644#issuecomment-597239761
!

module types

  implicit none
  private

  integer, parameter, public :: SP = selected_real_kind(6,30)
  integer, parameter, public :: DP = selected_real_kind(15,307)
  integer, parameter, public :: EP = selected_real_kind(18,90)
  integer, parameter, public :: QP = selected_real_kind(30,300)

  ! Working (FP) Precision
  integer, parameter, public :: WP = QP

end module types

module cbrt_lib
  use types, only: WP, SP, DP

  implicit none
  private

  ! Overloading
  interface cbrt
     module procedure cbrt_fike, cbrt_complex
  end interface

  public :: cbrt_apple, cbrt

contains

  ! Adapted from:
  ! https://people.freebsd.org/~lstewart/references/apple_tr_kt32_cuberoot.pdf
  function cbrt_apple(v) result (r)
    real(WP), intent(in) :: v
    real(WP) :: r

    real(WP), parameter :: C13 = 1.0_WP/3, C23 = 2.0_WP/3

    real(WP) :: fr, x
    integer :: ex, shx

    x = abs(v)

    ! Argument reduction
    !
    ! Separate into mantissa and exponent
    fr = fraction(x)
    ex = exponent(x)

    shx = mod(ex,3)

    ! Compute SHX such that (EX-SHX) is divisible by 3
    if (shx > 0) shx = shx-3

    ! Exponent of cube root
    ex = (ex-shx)/3

    ! 0.125 <= fr < 1.0
    ! SCALE(X,I) = X * RADIX(X)**I, the equivalent of C LDEXP(X,I)
    fr = scale(fr,shx)

    ! Compute seed with a quadratic approximation
    !
    ! 0.5 <= fr < 1
    fr = (-0.46946116_WP*fr+1.072302_WP)*fr+0.3812513_WP

    ! 6 bits of precision
    r = scale(fr,ex)

    ! Newton-Raphson iterations
    !
    ! 12 bits
    r = C23*r+C13*x/(r*r)

    ! 24 bits
    r = C23*r+C13*x/(r*r)

    if (WP == SP) then
       r = sign(r,v)
       return
    end if

    ! 48 bits
    r = C23*r+C13*x/(r*r)

    ! 96 bits
    r = C23*r+C13*x/(r*r)

    if (WP == DP) then
       r = sign(r,v)
       return
    end if

    ! 192 bits QP
    r = C23*r+C13*x/(r*r)

    r = sign(r,v)

  end function cbrt_apple

  ! Implemented from https://epdf.tips (Computer Evaluation of
  ! Mathematical Functions, C. T. FIKE) and suggestions extracted from
  ! above (apple_tr_kt32_cuberoot.pdf)
  !
  function cbrt_fike(v) result (y)
    real(WP), intent(in) :: v
    real(WP) :: y

    ! C16 = 1/(16)**(1/3), C256 = 1/(256)**(1/3)
    real(WP), parameter :: A = 0.391807_WP, B = 0.706799_WP, &
         C16 = 0.396850262992049869_WP, &
         C256 = 0.157490131236859146_WP

    real(WP) :: m, x
    integer :: q, r

    x = abs(v)

    ! Argument reduction
    !
    ! Separate into mantissa and exponent
    m = fraction(x)
    q = exponent(x)

    r = mod(q,4)

    ! Compute R such that (Q-R) is divisible by 4
    if (r > 0) r = r-4

    ! Exponent of cube root
    q = (q-r)/4

    ! 0.0625 <= m < 1.0  (1/16 == 0.625)
    ! SCALE(X,I) = X * RADIX(X)**I, the equivalent of C LDEXP(X,I)
    m = scale(m,r)

    ! Compute seed with a linear approximation (r(mp)). See FIKE book
    m = A+B*m

    r = mod(q,3)

    select case (r)
    case (0)
       ! 16 ** (q/3)r(mp)
       y = scale(m,4*(q/3))
    case (2)
       ! 16 ** ((q+1)/3)[r(mp)/16**(1/3)]
       y = scale(m*C16,4*((q+1)/3))
    case (1)
       ! 16 ** ((q+2)/3)[r(mp)/256**(1/3)]
       y = scale(m*C256,4*((q+2)/3))
    case default
       ! 16 ** (q/3)r(mp)
       y = scale(m,4*(q/3))
    end select

    ! Newton-Raphson iterations
    !
    ! RHO(0) = 2 ** (-3.34)
    ! RHO(k+1) = (2/3)*RHO(k) ** 3,  k = 0, 1, 2, ...:
    !
    !   RHO(1) = 2 ** (-10.605)
    !   RHO(2) = 2 ** (-32.400)
    !   RHO(3) = 2 ** (-97.785)
    !   RHO(4) = 2 ** (-293.940)
    !
    ! RHO is, say, the tolerance (DELTA X/ X)
    !

    ! First iteration (say 10 bit)
    m = y*y
    y = y-(m-(x/y))/(2*y+x/m)

    ! Second iteration (say 32 bit)
    m = y*y
    y = y-(m-(x/y))/(2*y+x/m)

    ! Third iteration (say 97 bit)
    m = y*y
    y = y-(m-(x/y))/(2*y+x/m)

    ! Fourth iteration (say 293 bit)
    m = y*y
    y = y-(m-(x/y))/(2*y+x/m)

    y = sign(y,v)

  end function cbrt_fike

  ! Return ALL the three complex roots
  function cbrt_complex(z) result(w)
    complex(WP), intent(in) :: z
    complex(WP) :: w(3)

    real(WP), parameter :: S3 = sqrt(3.0_WP)

    real(WP) :: phi, a, b, r, a3, b3

    a = real(z)
    b = aimag(z)

    ! arg(z)/3
    phi = atan2(b,a)/3

    ! CBRT(R)
    r = abs(z)
    r = cbrt_fike(r)

    a = r*cos(phi)
    b = r*sin(phi)

    w(1) = cmplx(a,b,kind=WP)

    a = -0.5_WP*a
    b = -0.5_WP*b

    a3 = S3*a
    b3 = S3*b

    w(2) = cmplx(a+b3,b-a3,kind=WP)
    w(3) = cmplx(a-b3,b+a3,kind=WP)

  end function cbrt_complex
end module cbrt_lib

program cbrt_test
  use types, only: WP
  use cbrt_lib, only: cbrt_apple, cbrt
  implicit none

  integer :: i
  real(WP) :: x, x3, fx3

  print *, cbrt_apple(1.0_WP), cbrt_apple(-1.0_WP)
  print *, cbrt_apple(8.0_WP), cbrt_apple(-8.0_WP)
  print *, cbrt_apple(27.0_WP), cbrt_apple(-27.0_WP)
  print *, cbrt_apple(64.0_WP), cbrt_apple(-64.0_WP)

  ! This app:
  !   4.32674871092222514696491493234032892 (QP)
  !
  ! Online (https://keisan.casio.com/calculator):
  !   4.3267487109222251469649149323403287652 (38 digits)
  print *, cbrt_apple(81.0_WP), cbrt_apple(-81.0_WP)
  print *, cbrt_apple(125.0_WP), cbrt_apple(-125.0_WP)
  print *, cbrt_apple(1000.0_WP), cbrt_apple(-1000.0_WP)

  print *
  print *, cbrt(1.0_WP), cbrt(-1.0_WP)
  print *, cbrt(8.0_WP), cbrt(-8.0_WP)
  print *, cbrt(27.0_WP), cbrt(-27.0_WP)
  print *, cbrt(64.0_WP), cbrt(-64.0_WP)
  print *, cbrt(81.0_WP), cbrt(-81.0_WP)
  print *, cbrt(125.0_WP), cbrt(-125.0_WP)
  print *, cbrt(1000.0_WP), cbrt(-1000.0_WP)

  ! Adapted from:
  !
  ! https://github.com/fortran-lang/stdlib/issues/214#issuecomment-656625976
  !
  do i = 1, 100
     call random_number(x)
     x = x*10.0_WP**i
     x3 = cbrt_apple(x)
     fx3 = cbrt(x)
     write(*,*) i, x, abs(x-x3**3)/x, abs(x-fx3**3)/x
  end do

  print *

  ! https://github.com/symengine/symengine/pull/1644
  print *, cbrt((-8.0_WP,0))

  ! CBRT() example 2 from IMSL, https://help.imsl.com/fortran/current/pdf/Fortran_Special_Functions_Library.pdf, page 19
  print *, cbrt((-3.0_WP,0.0076_WP))

  print *, cbrt((1.0_WP,0.0_WP))
  print *, cbrt((-1.0_WP,0.0_WP))
end program cbrt_test