A negative base raised to a rational exponent should not produce a co… by Geropy · Pull Request #1644 · symengine/symengine

Geropy · 2020-03-05T17:30:52Z

…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

…mplex number if the denominator of the exponent is odd; the result should be real. Added is_even() method to Integer to help accomplish this

isuruf · 2020-03-05T18:05:10Z

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 · 2020-03-05T18:21:51Z

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 · 2020-03-05T18:21:53Z

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 · 2020-03-05T18:22:07Z

closed by accident

isuruf · 2020-03-10T17:50:02Z

ping @certik

certik · 2020-03-10T18:21:41Z

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 · 2020-03-10T19:04:56Z

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 · 2020-03-10T19:27:38Z

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 · 2020-03-10T20:00:28Z

@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 · 2020-03-11T12:40:37Z

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

angelog0 · 2022-06-10T12:53:09Z

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

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

ca43741

…mplex number if the denominator of the exponent is odd; the result should be real. Added is_even() method to Integer to help accomplish this

Geropy closed this Mar 5, 2020

Geropy reopened this Mar 5, 2020

certik mentioned this pull request Jun 24, 2020

Cube root of a real number fortran-lang/stdlib#214

Open

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

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

A negative base raised to a rational exponent should not produce a co…#1644
Geropy wants to merge 1 commit into
symengine:masterfrom
Geropy:powreal-fix

Geropy commented Mar 5, 2020

Uh oh!

isuruf commented Mar 5, 2020

Uh oh!

Geropy commented Mar 5, 2020

Uh oh!

Geropy commented Mar 5, 2020

Uh oh!

Geropy commented Mar 5, 2020

Uh oh!

isuruf commented Mar 10, 2020

Uh oh!

certik commented Mar 10, 2020 •

edited

Loading

Uh oh!

certik commented Mar 10, 2020 •

edited

Loading

Uh oh!

Geropy commented Mar 10, 2020

Uh oh!

certik commented Mar 10, 2020

Uh oh!

Geropy commented Mar 11, 2020

Uh oh!

angelog0 commented Jun 10, 2022 •

edited

Loading

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

Conversation

Geropy commented Mar 5, 2020

Uh oh!

isuruf commented Mar 5, 2020

Uh oh!

Geropy commented Mar 5, 2020

Uh oh!

Geropy commented Mar 5, 2020

Uh oh!

Geropy commented Mar 5, 2020

Uh oh!

isuruf commented Mar 10, 2020

Uh oh!

certik commented Mar 10, 2020 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

certik commented Mar 10, 2020 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

Geropy commented Mar 10, 2020

Uh oh!

certik commented Mar 10, 2020

Uh oh!

Geropy commented Mar 11, 2020

Uh oh!

angelog0 commented Jun 10, 2022 • edited Loading Uh oh! There was an error while loading. Please reload this page.

Uh oh!

Uh oh!

Reviewers

Assignees

Labels

Projects

Milestone

Development

Uh oh!

4 participants

certik commented Mar 10, 2020 •

edited

Loading

certik commented Mar 10, 2020 •

edited

Loading

angelog0 commented Jun 10, 2022 •

edited

Loading