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

julia/base/complex.jl

Lines 875 to 878 in a0bb006

	function atan(z::Complex)
	w = atanh(Complex(-imag(z),real(z)))
	Complex(imag(w),-real(w))
	end

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

julia/complex.j

Line 220 in 44e0594

atanh(z::Complex) = log(sqrt((1+z)/(1-z)))

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
  

  julia/base/complex.jl

  Line 551 in a0bb006

  angle(z::Complex) = atan(imag(z), real(z))

• atan using atanh
  

  julia/base/complex.jl

  Lines 875 to 878 in a0bb006

  	function atan(z::Complex)
  	w = atanh(Complex(-imag(z),real(z)))
  	Complex(imag(w),-real(w))
  	end

• atanh using atan and angle
  

  julia/base/complex.jl

  Line 962 in a0bb006

  η = copysign(oftype(y,pi)/2 + atan(ym/2), y)/2

  julia/base/complex.jl

  Line 971 in a0bb006

  η = angle(Complex((1-x)*(1+x)-ysq, 2y))/2

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

julia/base/complex.jl

Line 955 in 795740a

elseif ax==1

is incorrect. If x == -1 and y != 0, then

julia/base/complex.jl

Lines 960 to 962 in 795740a

	ym = ay+ρ
	ξ = log(sqrt(sqrt(4+y*y))/sqrt(ym))
	η = copysign(oftype(y,pi)/2 + atan(ym/2), y)/2

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?