Julia: atan(1.0+im) is incorrectly signed

[JeffreySarnoff]

atan(1.0+im) == 1.0172219678978514 - 0.40235947810852507im

The sign of the imaginary part should be positive.

[inkydragon]

Note: tan(atan(z)) == z

julia> tan(atan(1+im))
1.0 - 1.0im

julia> tan(atan(1-im))
1.0 - 1.0im

julia> tan(atan(-1-im))
-1.0 - 1.0im

julia> tan(atan(-1+im))
-1.0 - 1.0im

julia> versioninfo()
Julia Version 1.0.3
Commit 099e826241 (2018-12-18 01:34 UTC)
Platform Info:
  OS: Windows (x86_64-w64-mingw32)
  CPU: Intel(R) Core(TM) i7-5600U CPU @ 2.60GHz
  WORD_SIZE: 64
  LIBM: libopenlibm
  LLVM: libLLVM-6.0.0 (ORCJIT, broadwell)

atan using atanh

https://github.com/JuliaLang/julia/blob/a0bb006d22945587598399c37b331003add4e912/base/complex.jl#L875-L878

we have

atan(z) == -im * atanh(im*z)

• Equ 4.4.22 @ Page 80 | Handbook of Mathematical Functions. by Abramowitz and Stegun.

z*im equal to complex(-z.img, z.real). So the upper code is right.

There must be something wrong with function atanh which has 40 lines code.

• atanh in julia/complex.jl · JuliaLang/julia

Also: tanh(atanh(z)) == z

julia> tanh(atanh(-1+im))
1.0 + 1.0im

julia> tanh(atanh(-1-im))
1.0 - 1.0im

julia> tanh(atanh(1-im))
1.0 - 1.0im

julia> tanh(atanh(1+im))
1.0 + 1.0im

---

A simple way arctan(z::Complex) = im/2.0 * ( log(im+z) - log(im-z) ) works fine.

julia> arctan(z::Complex) = im/2.0 * ( log(im+z) - log(im-z) )
arctan (generic function with 1 method)

julia> arctan(1+im)
1.0172219678978514 + 0.40235947810852507im

julia> tan(arctan(1+im))
1.0 + 1.0im

julia> tan(arctan(1-im))
1.0 - 1.0im

julia> tan(arctan(-1-im))
-1.0 - 1.0im

julia> tan(arctan(-1+im))
-1.0 + 1.0im

C# using this formula
ref:

• Equ 4.4.28 @ Page 80 | Handbook of Mathematical Functions. by Abramowitz and Stegun.
• corefx/Complex.cs · dotnet/corefx

[JeffreySarnoff]

julia> atan(1.0-im)
1.0172219678978514 - 0.40235947810852507im

julia> atan(1.0+im)
1.0172219678978514 - 0.40235947810852507im

julia> z = 1.0-im; tan(atan(z)) == z
true

julia> z = 1.0+im; tan(atan(z)) == z
false

With the sign of the imaginary part of atan(1.0±im) matching the sign of the imaginary part of 1.0±im, in both cases tan(atan(1.0±im)) == 1.0±im. This is the evaluation that Maple, Mathematica, Nemo use.

[chethega]

Quick test for some of the branch cuts:

julia> for r in [+1.0, -1.0] for op in [identity, nextfloat, prevfloat] for im_ in [0.0, nextfloat(0.0), prevfloat(0.0), 0.3]
       z = op(r) + im_*im
       @show z
       try
           @show tanh(atanh(z)), atanh(z)
       catch ex
       @show ex
       end end end end
z = 1.0 + 0.0im
(tanh(atanh(z)), atanh(z)) = (1.0 + 0.0im, Inf + 0.0im)
z = 1.0 + 5.0e-324im
(tanh(atanh(z)), atanh(z)) = (1.0 + 2.983336292480124e-154im, 177.09910463306602 + 0.7853981633974483im)
z = 1.0 - 5.0e-324im
(tanh(atanh(z)), atanh(z)) = (1.0 - 2.983336292480124e-154im, 177.09910463306602 - 0.7853981633974483im)
z = 1.0 + 0.3im
(tanh(atanh(z)), atanh(z)) = (1.0 + 0.3im, 0.9541226446766455 + 0.8598431372021968im)
z = 1.0000000000000002 + 0.0im
(tanh(atanh(z)), atanh(z)) = (1.0000000000000002 + 2.719262146893785e-32im, 18.36840028483855 + 1.5707963267948966im)
z = 1.0000000000000002 + 5.0e-324im
(tanh(atanh(z)), atanh(z)) = (1.0000000000000002 + 2.719262146893785e-32im, 18.36840028483855 + 1.5707963267948966im)
z = 1.0000000000000002 - 5.0e-324im
(tanh(atanh(z)), atanh(z)) = (1.0000000000000002 - 2.719262146893785e-32im, 18.36840028483855 - 1.5707963267948966im)
z = 1.0000000000000002 + 0.3im
(tanh(atanh(z)), atanh(z)) = (1.0000000000000004 + 0.30000000000000004im, 0.9541226446766456 + 0.8598431372021974im)
z = 0.9999999999999999 + 0.0im
(tanh(atanh(z)), atanh(z)) = (0.9999999999999998 + 0.0im, 18.714973875118524 + 0.0im)
z = 0.9999999999999999 + 5.0e-324im
(tanh(atanh(z)), atanh(z)) = (0.9999999999999998 + 5.0e-324im, 18.714973875118524 + 2.2250738585072014e-308im)
z = 0.9999999999999999 - 5.0e-324im
(tanh(atanh(z)), atanh(z)) = (0.9999999999999998 - 5.0e-324im, 18.714973875118524 - 2.2250738585072014e-308im)
z = 0.9999999999999999 + 0.3im
(tanh(atanh(z)), atanh(z)) = (1.0 + 0.3im, 0.9541226446766455 + 0.8598431372021968im)
z = -1.0 + 0.0im
(tanh(atanh(z)), atanh(z)) = (-1.0 + 0.0im, -Inf + 0.0im)
z = -1.0 + 5.0e-324im
(tanh(atanh(z)), atanh(z)) = (1.0 + 2.983336292480124e-154im, 177.09910463306602 + 0.7853981633974483im)
z = -1.0 - 5.0e-324im
(tanh(atanh(z)), atanh(z)) = (1.0 - 2.983336292480124e-154im, 177.09910463306602 - 0.7853981633974483im)
z = -1.0 + 0.3im
(tanh(atanh(z)), atanh(z)) = (1.0 + 0.3im, 0.9541226446766455 + 0.8598431372021968im) #wrong
z = -0.9999999999999999 + 0.0im
(tanh(atanh(z)), atanh(z)) = (-0.9999999789265759 + 0.0im, -9.184200142419275 + 0.0im)
z = -0.9999999999999999 + 5.0e-324im
(tanh(atanh(z)), atanh(z)) = (-0.9999999789265759 + 9.37798487e-316im, -9.184200142419275 + 2.2250738585072014e-308im)
z = -0.9999999999999999 - 5.0e-324im
(tanh(atanh(z)), atanh(z)) = (-0.9999999789265759 - 9.37798487e-316im, -9.184200142419275 - 2.2250738585072014e-308im)
z = -0.9999999999999999 + 0.3im
(tanh(atanh(z)), atanh(z)) = (-0.9999999999999999 + 0.3000000000000003im, -0.9541226446766451 + 0.8598431372021968im)
z = -1.0000000000000002 + 0.0im
ex = DomainError(-1.0000000000000002, "log1p will only return a complex result if called with a complex argument. Try log1p(Complex(x)).") #looks super wrong
z = -1.0000000000000002 + 5.0e-324im
ex = DomainError(-1.0000000000000002, "log1p will only return a complex result if called with a complex argument. Try log1p(Complex(x)).") #looks super wrong
z = -1.0000000000000002 - 5.0e-324im
ex = DomainError(-1.0000000000000002, "log1p will only return a complex result if called with a complex argument. Try log1p(Complex(x)).") #looks super wrong
z = -1.0000000000000002 + 0.3im
(tanh(atanh(z)), atanh(z)) = (-1.0000000000000007 + 0.29999999999999805im, -0.9541226446766489 + 0.8598431372021974im)

We maybe should have test-suites like this for more of the complex functions, possibly also including prevfloat(typemax(T)) and similar.

[inkydragon]

some history:

this version (oldest one, 9 year ago) is right
https://github.com/JuliaLang/julia/blob/44e0594da8029cf8acc49b29dfaa0fe50f53a07d/complex.j#L220

After pr #2891, next version is wrong
@jiahao @JeffBezanson

• Implements Kahan's versions of elementary complex functions… · JuliaLang/julia@60d266a

atanh_60d266

function atanh_60d266(z::Complex{T}) where T<:AbstractFloat
    Ω=prevfloat(typemax(T))
    θ=sqrt(Ω)/4
    ρ=1/θ
    x=real(z)
    y=imag(z)
    if x > θ || abs(y) > θ #Prevent overflow
        return complex(copysign(pi/2, y), real(1/z))
    elseif x==one(x)
        ym=abs(y)+ρ
        ξ=log(sqrt(sqrt(4+y^2))/sqrt(ym))
        η=copysign(pi/2+atan(ym/2), y)/2
    else #Normal case
        ysq=(abs(y)+ρ)^2
        ξ=log1p(4x/((1-x)^2 + ysq))/4
        η=angle(complex(((1-x)*(1+x)-ysq)/2, y))
    end
    complex(ξ, η)
end

julia> tanh(atanh_60d266(-1.0+im))
-1.2060113295832982 - 1.078689325833263im

julia> tanh(atanh_60d266(-1.0-im))
-1.2060113295832982 + 1.078689325833263im

and then all those later versions are wrong.

pr #2891 have mentioned a paper: Branch Cuts for Complex Elementary Functions (or: Much Ado About Nothing's Sign Bit)

We may need to take a closer look at this paper to fix this bug.

[JeffreySarnoff]

good digging -- I remember that we had it right at some point. That is the seminal paper.

[JeffreySarnoff]

Here is some function specific information on branch cuts generated from Maple.
Complex arctrig(h) functions according to Maple
Principal branches according to Maple

[inkydragon]

Some debug notes:

• angle using atan
  https://github.com/JuliaLang/julia/blob/a0bb006d22945587598399c37b331003add4e912/base/complex.jl#L551
• atan using atanh
  https://github.com/JuliaLang/julia/blob/a0bb006d22945587598399c37b331003add4e912/base/complex.jl#L875-L878
• atanh using atan and angle
  https://github.com/JuliaLang/julia/blob/a0bb006d22945587598399c37b331003add4e912/base/complex.jl#L962 https://github.com/JuliaLang/julia/blob/a0bb006d22945587598399c37b331003add4e912/base/complex.jl#L971

if you want to test atanh, you may choose a Correct definition of atan and angle first.

[ViralBShah]

We definitely need to get this right.

[JeffreySarnoff]

Wolfram Research has this well curated resource.

use the show all button on the upper left after selecting a function specific subtopic for example, this

I have added a pdf with Maple's plots of the principal branches here

[cafaxo]

It seems like
https://github.com/JuliaLang/julia/blob/795740a3948649f74fe1070c870c0afbb4501bec/base/complex.jl#L955
is incorrect. If x == -1 and y != 0, then
https://github.com/JuliaLang/julia/blob/795740a3948649f74fe1070c870c0afbb4501bec/base/complex.jl#L960-L962
is executed, but that code is only valid for x == 1.
(Reference: https://people.freebsd.org/~das/kahan86branch.pdf)

[JeffreySarnoff]

We definitely need to get this right.

I found this state of affairs obliquely. It is necessary to re-vet all the complex arc functions.

[simonbyrne]

Thanks for looking into it @cafaxo. Care to open a pull request?

[inkydragon]

Patch

change those judgment statement did fix the sign bug.
@cafaxo

Sign Test Result

function test_atanh()
    float_data_list = [-1.0, 0.0, 1.0]
    op_list = [identity, nextfloat, prevfloat]

    for op in op_list
    for real in float_data_list
    for imag in float_data_list
        z = op(real) + op(imag)*im
        @show z
        try 
            if tanh(atanh(z)) == z
                nothing
            elseif tanh(atanh(z)) ≈ z
                nothing
            else
                @show tanh(atanh(z)) atanh(z)
            end
        catch ex
            @show ex
        end
    end 
    end 
    end

end

julia> test_f()
z = -1.0 - 1.0im
z = -1.0 + 0.0im
z = -1.0 + 1.0im
z = 0.0 - 1.0im
z = 0.0 + 0.0im
z = 0.0 + 1.0im
z = 1.0 - 1.0im
z = 1.0 + 0.0im
z = 1.0 + 1.0im
z = -0.9999999999999999 - 0.9999999999999999im
z = -0.9999999999999999 + 5.0e-324im
tanh(atanh(z)) = -0.9999999789265759 + 9.37798487e-316im
atanh(z) = -9.184200142419275 + 2.2250738585072014e-308im
z = -0.9999999999999999 + 1.0000000000000002im
z = 5.0e-324 - 0.9999999999999999im
z = 5.0e-324 + 5.0e-324im
z = 5.0e-324 + 1.0000000000000002im
z = 1.0000000000000002 - 0.9999999999999999im
z = 1.0000000000000002 + 5.0e-324im
z = 1.0000000000000002 + 1.0000000000000002im
z = -1.0000000000000002 - 1.0000000000000002im
z = -1.0000000000000002 - 5.0e-324im
ex = DomainError(-1.0000000000000002, "log1p will only return a complex result if called with a complex argument. Try log1p(Complex(x)).")
z = -1.0000000000000002 + 0.9999999999999999im
z = -5.0e-324 - 1.0000000000000002im
z = -5.0e-324 - 5.0e-324im
z = -5.0e-324 + 0.9999999999999999im
z = 0.9999999999999999 - 1.0000000000000002im
z = 0.9999999999999999 - 5.0e-324im
z = 0.9999999999999999 + 0.9999999999999999im

BUT it cause test about NaN, Inf failed.

 base/complex.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/base/complex.jl b/base/complex.jl
index dc736ebb..6f0006dc 100644
--- a/base/complex.jl
+++ b/base/complex.jl
@@ -905,7 +905,7 @@ function atanh(z::Complex{T}) where T<:AbstractFloat
     x, y = reim(z)
     ax = abs(x)
     ay = abs(y)
-    if ax > θ || ay > θ #Prevent overflow
+    if x > θ || ay > θ #Prevent overflow
         if isnan(y)
             if isinf(x)
                 return Complex(copysign(zero(x),x), y)
@@ -917,7 +917,7 @@ function atanh(z::Complex{T}) where T<:AbstractFloat
             return Complex(copysign(zero(x),x), copysign(oftype(y,pi)/2, y))
         end
         return Complex(real(1/z), copysign(oftype(y,pi)/2, y))
-    elseif ax==1
+    elseif x == 1
         if y == 0
             ξ = copysign(oftype(x,Inf),x)
             η = zero(y)

After apply this patch, all those if-else statement are the same as in the paper.
Although return statement still has a little difference：

Complex(ξ, η)
copysign(1, x) * Complex(ξ, -1*η)

Test

include("..\test\complex.jl") Result

Test Summary:                       |     Pass  Fail  Error  Broken     Total
  Overall                           | 37529466    13      1  327825  37857305
    ...
    complex                         |     8238    12              2      8252
    ...

julia> using Test

julia> include("..\\test\\complex.jl")
atanh(op(z)) == op(atanh(z)) for op in (conj, -): Test Failed at ~\julia-1.0.3\test\complex.jl:699
  Expression: isequal(atanh(complex(-1.0, 0.0)), complex(-Inf, 0.0))
   Evaluated: isequal(-Inf + 1.5707963267948966im, -Inf + 0.0im)
atanh(op(z)) == op(atanh(z)) for op in (conj, -): Test Failed at ~\julia-1.0.3\test\complex.jl:715
  Expression: isequal(atanh(complex(-Inf, 0.0)), complex(-0.0, pi / 2))
   Evaluated: isequal(NaN + 1.5707963267948966im, -0.0 + 1.5707963267948966im)
atanh(op(z)) == op(atanh(z)) for op in (conj, -): Test Failed at ~\julia-1.0.3\test\complex.jl:716
  Expression: isequal(atanh(complex(-Inf, -0.0)), complex(-0.0, -pi / 2))
   Evaluated: isequal(NaN - 1.5707963267948966im, -0.0 - 1.5707963267948966im)
atanh(op(z)) == op(atanh(z)) for op in (conj, -): Test Failed at ~\julia-1.0.3\test\complex.jl:717
  Expression: isequal(atanh(complex(-Inf, 5.0)), complex(-0.0, pi / 2))
   Evaluated: isequal(NaN + 1.5707963267948966im, -0.0 + 1.5707963267948966im)
atanh(op(z)) == op(atanh(z)) for op in (conj, -): Test Failed at ~\julia-1.0.3\test\complex.jl:718
  Expression: isequal(atanh(complex(-Inf, -5.0)), complex(-0.0, -pi / 2))
   Evaluated: isequal(NaN - 1.5707963267948966im, -0.0 - 1.5707963267948966im)
atanh(op(z)) == op(atanh(z)) for op in (conj, -): Test Failed at ~\julia-1.0.3\test\complex.jl:721
  Expression: isequal(atanh(complex(-Inf, NaN)), complex(-0.0, NaN))
   Evaluated: isequal(NaN + NaN*im, -0.0 + NaN*im)
atan(z) == -i*atanh(iz): Test Failed at ~\julia-1.0.3\test\complex.jl:735
  Expression: isequal(atan(complex(0.0, 1.0)), complex(0.0, Inf))
   Evaluated: isequal(1.5707963267948966 + Inf*im, 0.0 + Inf*im)
atan(z) == -i*atanh(iz): Test Failed at ~\julia-1.0.3\test\complex.jl:736
  Expression: isequal(atan(complex(0.0, Inf)), complex(pi / 2, 0.0))
   Evaluated: isequal(1.5707963267948966 + NaN*im, 1.5707963267948966 + 0.0im)
atan(z) == -i*atanh(iz): Test Failed at ~\julia-1.0.3\test\complex.jl:742
  Expression: isequal(atan(complex(-0.0, Inf)), complex(-pi / 2, 0.0))
   Evaluated: isequal(-1.5707963267948966 + NaN*im, -1.5707963267948966 + 0.0im)
atan(z) == -i*atanh(iz): Test Failed at ~\julia-1.0.3\test\complex.jl:746
  Expression: isequal(atan(complex(5.0, Inf)), complex(pi / 2, 0.0))
   Evaluated: isequal(1.5707963267948966 + NaN*im, 1.5707963267948966 + 0.0im)
atan(z) == -i*atanh(iz): Test Failed at ~\julia-1.0.3\test\complex.jl:750
  Expression: isequal(atan(complex(-5.0, Inf)), complex(-pi / 2, 0.0))
   Evaluated: isequal(-1.5707963267948966 + NaN*im, -1.5707963267948966 + 0.0im)
atan(z) == -i*atanh(iz): Test Failed at ~\julia-1.0.3\test\complex.jl:774
  Expression: isequal(atan(complex(NaN, Inf)), complex(NaN, 0.0))
   Evaluated: isequal(NaN + NaN*im, NaN + 0.0im)
Test.DefaultTestSet("more cpow", Any[Test Broken
  Expression: (Inf + 1im) ^ 3 === (Inf + 1im) ^ 3.0 === (Inf + 1im) ^ (3 + 0im) === Inf + Inf * im, Test Broken
  Expression: (Inf + 1im) ^ 3.1 === (Inf + 1im) ^ (3.1 + 0im) === Inf + Inf * im], 7125, false)

TODO

• [ ] find out why there is a DomainError about log1p
• [ ] fix other error in test_atanh()
• [ ] fix failed NaN and Inf test
• [ ] add more test case to \test\complex.jl

[JeffreySarnoff]

is there a ready explanation of why the returned expressions differ
Complex(ξ, η)
copysign(1, x) * Complex(ξ, -1*η)

[JeffreySarnoff]

This 2017 paper gives values used for testing the implementation of complex elementary functions
On quality of implementation of Fortran 2008 complex intrinsic functions on branch cuts

[cafaxo]

I suggest

function atanh(z::Complex{T}) where T<:AbstractFloat
    Ω = prevfloat(typemax(T))
    θ = sqrt(Ω)/4
    ρ = 1/θ
    x, y = reim(z)
    ax = abs(x)
    ay = abs(y)
    if ax > θ || ay > θ #Prevent overflow
        if isnan(y)
            if isinf(x)
                return Complex(copysign(zero(x),x), y)
            else
                return Complex(real(1/z), y)
            end
        end
        if isinf(y)
            return Complex(copysign(zero(x),x), copysign(oftype(y,pi)/2, y))
        end
        return Complex(real(1/z), copysign(oftype(y,pi)/2, y))
    end

    β = copysign(1, x)
    z = β * conj(z)
    x, y = reim(z)

    if x == 1
        if y == 0
            ξ = oftype(x, Inf)
            η = y
        else
            ym = ay+ρ
            ξ = log(sqrt(sqrt(4+y*y))/sqrt(ym))
            η = copysign(oftype(y,pi)/2 + atan(ym/2), y)/2
        end
    else #Normal case
        ysq = (ay+ρ)^2
        if x == 0
            ξ = x
        else
            ξ = log1p(4x/((1-x)^2 + ysq))/4
        end
        η = angle(Complex((1-x)*(1+x)-ysq, 2y))/2
    end

    return β * complex(ξ, -η)
end

As a patch

@@ -952,10 +952,16 @@
             return Complex(copysign(zero(x),x), copysign(oftype(y,pi)/2, y))
         end
         return Complex(real(1/z), copysign(oftype(y,pi)/2, y))
-    elseif ax==1
+    end
+
+    β = copysign(1, x)
+    z = β * conj(z)
+    x, y = reim(z)
+
+    if x == 1
         if y == 0
-            ξ = copysign(oftype(x,Inf),x)
-            η = zero(y)
+            ξ = oftype(x, Inf)
+            η = y
         else
             ym = ay+ρ
             ξ = log(sqrt(sqrt(4+y*y))/sqrt(ym))
@@ -970,7 +976,8 @@
         end
         η = angle(Complex((1-x)*(1+x)-ysq, 2y))/2
     end
-    Complex(ξ, η)
+    
+    return β * complex(ξ, -η)
 end
 atanh(z::Complex) = atanh(float(z))

This fixes the sign issue, passes all tests in Base, passes test_atanh from @inkydragon, and does not throw on atanh(prevfloat(-1.0) + 0im). It also matches the paper more closely.

[JeffreySarnoff]

@cafaxo thank you for the diligence and attention

[cafaxo]

I do not know why the paper conjugates z. This should probably be dropped from my suggestion.

[JeffreySarnoff]

hold off on that
conj(complex) has IEEE754 specific details as I recall
check what conj does with re aor im parts +0, -0

Jeffrey Sarnoff

On Feb 14, 2019, at 10:57 AM, cafaxo [e-mail] wrote:

I do not know why the paper conjugates z. This should probably be dropped from my suggestion.

—
You are receiving this because you were mentioned.
Reply to this email directly, view it on GitHub, or mute the thread.

[cafaxo]

conj only negates the imaginary part:
https://github.com/JuliaLang/julia/blob/bd07ef4f97f3ab00d56719eb7c38bfce2ea43344/base/complex.jl#L259
I am pretty sure that we do not need that here.

[JeffreySarnoff]

exerpted

x, y = reim(z) # ... β = copysign(1, x) z = β * conj(z) x, y = reim(z)
z = β * conj(z) is not z = β * z

julia> z = Complex(1.0, +0.0); copysign(1, real(z)) * conj(z)
1.0 - 0.0im

julia> z = Complex(1.0, +0.0); copysign(1, real(z)) * (z)
1.0 + 0.0im

julia> z = Complex(1.0, -0.0); copysign(1, real(z)) * conj(z)
1.0 + 0.0im

julia> z = Complex(1.0, -0.0); copysign(1, real(z)) * (z)
1.0 - 0.0im

[simonbyrne]

Fixed by #31061