symengine
A negative base raised to a rational exponent should not produce a co…
…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
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?
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
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
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.
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.
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
@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?
In the "Possible issues" of the Surd definition, it says that it is undefined for even n's and negative real x
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
