Fix incorrect sign of atanh(complex(x,y)) if x == -1 #31061

cafaxo · 2019-02-13T19:01:15Z

In the case that x==-1, we have to flip the sign of ξ. Fixes the problem noticed in #31054.

Update:
The current Julia implementation was mistranscribed from the reference: The multiplication with the sign of x is missing (see #31061 (comment)). This PR currently adds that multiplication, fixing both the original sign issue and the problems (DomainError) near z=-1+0im.

In the case that x==-1, we have to flip the sign of ξ.

base/complex.jl

Co-Authored-By: cafaxo <cafaxo@gmail.com>

stevengj · 2019-02-13T22:11:57Z

The η = zero(y) in the if y == 0 also seems wrong.

Keno · 2019-02-13T23:00:58Z

This should have a test to prevent this from regressing.

StefanKarpinski · 2019-02-13T23:33:13Z

The η = zero(y) in the if y == 0 also seems wrong.

Should it be η = y in that case?

chethega · 2019-02-14T12:18:39Z

Can we also fix atanh(prevfloat(-1.0) + 0im)?

Most complex functions return some complex mix of Inf and NaN for singularities, instead of throwing. Worst case, NaN + NaN*im seems like a better alternative.

martinholters · 2019-02-14T13:07:32Z

I think we have at some point decided that throwing an appropriate error is more Julian than silently producing NaN. The exception is where something built-in gives NaN and, for performance, we don't do an extra check to throw instead. (0.0/0.0 would be the prime example for the latter.)

stevengj · 2019-02-14T14:38:12Z

Arguably in the y==0 case it should be η = copysign(oftype(y,pi)/4, y), but I guess we wanted to agree with atan(real(z)) in that case? But at least we should do η = y to not discard the sign of zero.

stevengj · 2019-02-14T14:41:10Z

Note that

julia> ccall(:catanh, ComplexF64, (ComplexF64,), complex(1.0, +0.0))
Inf + 0.7853981633974483im

julia> ccall(:catanh, ComplexF64, (ComplexF64,), complex(1.0, -0.0))
Inf - 0.7853981633974483im

so that would also argue for η = copysign(oftype(y,pi)/4, y) in the y==0 case.

(Python cmath.atanh throws a domain error here.)

stevengj · 2019-02-14T14:46:31Z

Note that this code was from #2891 by @jiahao, which cites Kahan, but I think it mistranscribed the reference, which says:

stevengj · 2019-02-14T14:49:10Z

Note that this reference doesn't have the copysign on the ξ term, but it multiplies by the sign at the end via β … it went through some extra machinations to deal with unsigned zeros on ancient hardware, apparently?

cafaxo · 2019-02-14T14:53:00Z

Yes, is differs in another way: It does z = β * conj(z). Thus it only has to handle the case x == 1. The case x == -1 cannot occur, in contrast to the current julia implementation.

cafaxo · 2019-02-14T17:28:22Z

We can fix atanh(prevfloat(-1.0) + 0im) by multiplying z with the sign of x and then multiplying (ξ, η) with that sign, as in the paper. As we do not have to handle x == -1 anymore, we can drop the copysign that this PR would currently introduce to fix the sign issue. This results in

julia> ccall(:catanh, ComplexF64, (ComplexF64,), prevfloat(-1.0) + 0im)
-18.36840028483855 + 1.5707963267948966im

julia> atanh(prevfloat(-1.0) + 0im)
-18.36840028483855 + 1.5707963267948966im

Should I add this change?

test/complex.jl

stevengj · 2019-02-14T19:02:25Z

I don't understand the rationale for the ay+ρ logic in the code (either Julia's or Kahan's). Eliminating that seems to make the code more accurate as well as simpler.

In particular, consider the following example where our current code gives quite inaccurate results:

julia> atanh(1+2.98e-154im)
176.7528106982835 + 0.7853981633974483im

julia> tanh(atanh(1+2.98e-154im))
1.0 + 5.963336292480018e-154im

In contrast, catanh(z) = ccall(:catanh, ComplexF64, (ComplexF64,), z) gives accurate results

julia> catanh(1+2.98e-154im)
177.09966410056285 + 0.7853981633974483im

julia> tanh(catanh(1+2.98e-154im))
1.0 + 2.9799999999999177e-154im

And if I simplify our atanh implementation to just eliminate the +ρ factors, I get similarly accurate results:

function myatanh(z::Complex{T}) where T<:AbstractFloat
    Ω = prevfloat(typemax(T))
    θ = sqrt(Ω)/4
    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))
    elseif ax==1
        ξ = copysign(log(sqrt(sqrt(4+y*y))/sqrt(ay)), x)
        η = copysign(oftype(y,pi)/2 + atan(ay/2), y)/2
    else #Normal case
        ysq = ay^2
        ξ = log1p(4x/((1-x)^2 + ysq))/4
        η = angle(Complex((1-x)*(1+x)-ysq, 2y))/2
    end
    Complex(ξ, η)
end

julia> myatanh(1+2.98e-154im)
177.09966410056285 + 0.7853981633974483im

julia> tanh(myatanh(1+2.98e-154im))
1.0 + 2.9799999999999177e-154im

What problem does the +ρ supposedly solve?

stevengj · 2019-02-14T19:07:36Z

As a more minor matter, I also don't understand why we do log(sqrt(sqrt(4+y*y)) rather than log(4+y*y)/4

stevengj · 2019-02-14T19:13:49Z

@cafaxo, I agree that flipping the sign seems like a good idea.

StefanKarpinski · 2019-02-14T20:56:32Z

This seems like a minor change to everyone on triage: technically breaking, but better, so we should change it but not backport that change. Keep in mind that users of long term support versions are super conservative and want their programs to keep working as is unless the previous behavior was so horribly wrong that it couldn't be working right.

JeffBezanson · 2019-02-15T02:43:58Z

Ok, triage might be wrong about this one.

chethega · 2019-02-15T21:34:01Z

Can't we just use a different expansion near -1.0+0im? if abs(x+1)+abs(y) < 1e-2 ξ = log((1+x)^2 +ysq)/((1-x)^2 +ysq)/4 end, or something.

Our expansion near that singularity is needlessly bad, which leads to rounding causing the DomainError. log1p is definitely the right expansion near z=0, and looks ok near z=1+0im.

cafaxo · 2019-02-15T22:16:25Z

The current Julia implementation mistranscribed the reference: The multiplication with the sign of x is missing (see #31061 (comment)). This PR currently adds that multiplication, fixing both the original sign issue and the problems near z=-1+0im.

StefanKarpinski · 2019-02-16T20:46:00Z

Should this also add some tests of for x ± 0.0im values to make sure that the branch cuts are correct?

simonbyrne · 2019-02-17T00:57:48Z

There are a bunch of branch cut and limit tests here:

julia/test/complex.jl

Lines 684 to 730 in dfee945

    
           @testset "atanh(op(z)) == op(atanh(z)) for op in (conj, -)" begin 
        
               @test isequal(atanh(complex( 0, 0)),complex( 0.0, 0.0)) #integer fallback 
        
               @test isequal(atanh(complex( 0.0, 0.0)),complex( 0.0, 0.0)) 
        
               @test isequal(atanh(complex( 0.0,-0.0)),complex( 0.0,-0.0)) 
        
               @test isequal(atanh(complex( 0.0, NaN)),complex( 0.0, NaN)) 
        
               @test isequal(atanh(complex( 0.0, Inf)),complex( 0.0, pi/2)) 
        
               @test isequal(atanh(complex( 0.0,-Inf)),complex( 0.0,-pi/2)) 
        
               @test isequal(atanh(complex(-0.0, NaN)),complex(-0.0, NaN)) 
        
               @test isequal(atanh(complex(-0.0, 0.0)),complex(-0.0, 0.0)) 
        
               @test isequal(atanh(complex(-0.0,-0.0)),complex(-0.0,-0.0)) 
        
               @test isequal(atanh(complex(-0.0, Inf)),complex(-0.0, pi/2)) 
        
               @test isequal(atanh(complex(-0.0,-Inf)),complex(-0.0,-pi/2)) 
        
               @test isequal(atanh(complex( 1.0, 0.0)),complex( Inf, 0.0)) 
        
               @test isequal(atanh(complex(-1.0, 0.0)),complex(-Inf, 0.0)) 
        
               @test isequal(atanh(complex( 5.0, Inf)),complex( 0.0, pi/2)) 
        
               @test isequal(atanh(complex( 5.0,-Inf)),complex( 0.0,-pi/2)) 
        
               @test isequal(atanh(complex( 5.0, NaN)),complex( NaN, NaN)) 
        
               @test isequal(atanh(complex(-5.0, Inf)),complex(-0.0, pi/2)) 
        
               @test isequal(atanh(complex(-5.0,-Inf)),complex(-0.0,-pi/2)) 
        
               @test isequal(atanh(complex(-5.0, NaN)),complex( NaN, NaN)) 
        
               @test isequal(atanh(complex( Inf, 0.0)),complex(0.0, pi/2)) 
        
               @test isequal(atanh(complex( Inf,-0.0)),complex(0.0,-pi/2)) 
        
               @test isequal(atanh(complex( Inf, 5.0)),complex(0.0, pi/2)) 
        
               @test isequal(atanh(complex( Inf,-5.0)),complex(0.0,-pi/2)) 
        
               @test isequal(atanh(complex( Inf, Inf)),complex(0.0, pi/2)) 
        
               @test isequal(atanh(complex( Inf,-Inf)),complex(0.0,-pi/2)) 
        
               @test isequal(atanh(complex( Inf, NaN)),complex(0.0, NaN)) 
        
               @test isequal(atanh(complex(-Inf, 0.0)),complex(-0.0, pi/2)) 
        
               @test isequal(atanh(complex(-Inf,-0.0)),complex(-0.0,-pi/2)) 
        
               @test isequal(atanh(complex(-Inf, 5.0)),complex(-0.0, pi/2)) 
        
               @test isequal(atanh(complex(-Inf,-5.0)),complex(-0.0,-pi/2)) 
        
               @test isequal(atanh(complex(-Inf, Inf)),complex(-0.0, pi/2)) 
        
               @test isequal(atanh(complex(-Inf,-Inf)),complex(-0.0,-pi/2)) 
        
               @test isequal(atanh(complex(-Inf, NaN)),complex(-0.0, NaN)) 
        
               @test isequal(atanh(complex( NaN, 0.0)),complex( NaN, NaN)) 
        
               @test isequal(atanh(complex( NaN,-0.0)),complex( NaN, NaN)) 
        
               @test isequal(atanh(complex( NaN, 5.0)),complex( NaN, NaN)) 
        
               @test isequal(atanh(complex( NaN,-5.0)),complex( NaN, NaN)) 
        
               @test isequal(atanh(complex( NaN, Inf)),complex( 0.0, pi/2)) 
        
               @test isequal(atanh(complex( NaN,-Inf)),complex( 0.0,-pi/2)) 
        
               @test isequal(atanh(complex( NaN, NaN)),complex( NaN, NaN)) 
        
           end

so can add those there.

StefanKarpinski · 2019-02-17T14:14:59Z

I guess they were not incorrect previously then. Carry on.

base/complex.jl

cafaxo · 2019-02-18T13:11:29Z

I have added tests for the cases that were fixed after replacing η = zero(y) with η = y.

simonbyrne · 2019-02-18T21:24:07Z

Looks like the CI failures are due to getting slightly different values here:

julia/stdlib/LinearAlgebra/src/dense.jl

Lines 927 to 930 in e5e1253

    
           julia> acos(cos([0.5 0.1; -0.2 0.3])) 
        
           2×2 Array{Complex{Float64},2}: 
        
             0.5-5.55112e-17im  0.1-2.77556e-17im 
        
            -0.2+2.498e-16im    0.3-3.46945e-16im

These should be updated with the new values.

simonbyrne · 2019-02-19T22:13:01Z

Thanks @cafaxo!

StefanKarpinski · 2019-02-20T19:20:19Z

It would be nice to clean up the commit message when squashing this kind of PR. The GitHub UI makes it quite easy to do—you just need to edit the message in the web interface.

simonbyrne · 2019-02-20T19:21:55Z

Good point, will do in future.

fp4code · 2019-03-27T13:52:23Z

A partial backport should be done! This is not about the choice of the cut, just
tanh(atanh(-1.0 + 0.1im)) ≈ 1.0 + 0.1im is a wrong result!

* Fix incorrect sign in atanh In the case that x==-1, we have to flip the sign of ξ. * Formatting: Add space after comma Co-Authored-By: cafaxo <cafaxo@gmail.com> * Add test * Do not drop sign of zero * Flip sign to avoid a DomainError This fixes atanh(prevfloat(-1.0) + 0im) * Test all four combinations Co-Authored-By: cafaxo <cafaxo@gmail.com> * Tests for expressions that were signed incorrectly * Use the correct type for 1 * Update doc * Update doc: Remove trailing whitespace (cherry picked from commit 66341a3)

Fix incorrect sign in atanh

c445183

In the case that x==-1, we have to flip the sign of ξ.

cafaxo changed the title ~~Fix incorrect sign in atanh if x == -1~~ Fix incorrect sign of atanh(complex(x,y)) if x == -1 Feb 13, 2019

StefanKarpinski reviewed Feb 13, 2019

View reviewed changes

base/complex.jl Outdated Show resolved Hide resolved

Formatting: Add space after comma

d987fed

Co-Authored-By: cafaxo <cafaxo@gmail.com>

fredrikekre added domain:maths Mathematical functions needs tests Unit tests are required for this change kind:bugfix This change fixes an existing bug status:triage This should be discussed on a triage call backport 1.0 labels Feb 13, 2019

Add test

b30bf57

fredrikekre removed the needs tests Unit tests are required for this change label Feb 13, 2019

simonbyrne reviewed Feb 14, 2019

View reviewed changes

test/complex.jl Outdated Show resolved Hide resolved

JeffBezanson removed backport 1.0 status:triage This should be discussed on a triage call labels Feb 14, 2019

stevengj reviewed Feb 18, 2019

View reviewed changes

base/complex.jl Outdated Show resolved Hide resolved

cafaxo added 2 commits February 18, 2019 14:00

Tests for expressions that were signed incorrectly

e3598fe

Use the correct type for 1

f8eacdb

cafaxo added 2 commits February 18, 2019 22:58

Update doc

f9cef28

Update doc: Remove trailing whitespace

e4f0326

stevengj approved these changes Feb 18, 2019

View reviewed changes

simonbyrne closed this Feb 18, 2019

simonbyrne reopened this Feb 18, 2019

simonbyrne merged commit 66341a3 into JuliaLang:master Feb 19, 2019

simonbyrne mentioned this pull request Feb 19, 2019

atan(1.0+im) is incorrectly signed #31054

Closed

simonbyrne added status:triage This should be discussed on a triage call backport 1.0 labels Feb 19, 2019

KristofferC removed the status:triage This should be discussed on a triage call label Aug 26, 2019

KristofferC mentioned this pull request Dec 3, 2019

WIP: Backports release 1.0.6 #34011

Closed

56 tasks

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

Fix incorrect sign of atanh(complex(x,y)) if x == -1 #31061

Fix incorrect sign of atanh(complex(x,y)) if x == -1 #31061

cafaxo commented Feb 13, 2019 •

edited

Loading

stevengj commented Feb 13, 2019

Keno commented Feb 13, 2019

StefanKarpinski commented Feb 13, 2019

chethega commented Feb 14, 2019

martinholters commented Feb 14, 2019

stevengj commented Feb 14, 2019 •

edited

Loading

stevengj commented Feb 14, 2019 •

edited

Loading

stevengj commented Feb 14, 2019

stevengj commented Feb 14, 2019 •

edited

Loading

cafaxo commented Feb 14, 2019

cafaxo commented Feb 14, 2019 •

edited

Loading

stevengj commented Feb 14, 2019 •

edited

Loading

stevengj commented Feb 14, 2019

stevengj commented Feb 14, 2019

StefanKarpinski commented Feb 14, 2019

JeffBezanson commented Feb 15, 2019

chethega commented Feb 15, 2019

cafaxo commented Feb 15, 2019 •

edited

Loading

StefanKarpinski commented Feb 16, 2019

simonbyrne commented Feb 17, 2019

StefanKarpinski commented Feb 17, 2019

cafaxo commented Feb 18, 2019 •

edited

Loading

simonbyrne commented Feb 18, 2019 •

edited

Loading

simonbyrne commented Feb 19, 2019

StefanKarpinski commented Feb 20, 2019

simonbyrne commented Feb 20, 2019

fp4code commented Mar 27, 2019 •

edited

Loading

Fix incorrect sign of atanh(complex(x,y)) if x == -1 #31061

Fix incorrect sign of atanh(complex(x,y)) if x == -1 #31061

Conversation

cafaxo commented Feb 13, 2019 • edited Loading

stevengj commented Feb 13, 2019

Keno commented Feb 13, 2019

StefanKarpinski commented Feb 13, 2019

chethega commented Feb 14, 2019

martinholters commented Feb 14, 2019

stevengj commented Feb 14, 2019 • edited Loading

stevengj commented Feb 14, 2019 • edited Loading

stevengj commented Feb 14, 2019

stevengj commented Feb 14, 2019 • edited Loading

cafaxo commented Feb 14, 2019

cafaxo commented Feb 14, 2019 • edited Loading

stevengj commented Feb 14, 2019 • edited Loading

stevengj commented Feb 14, 2019

stevengj commented Feb 14, 2019

StefanKarpinski commented Feb 14, 2019

JeffBezanson commented Feb 15, 2019

chethega commented Feb 15, 2019

cafaxo commented Feb 15, 2019 • edited Loading

StefanKarpinski commented Feb 16, 2019

simonbyrne commented Feb 17, 2019

StefanKarpinski commented Feb 17, 2019

cafaxo commented Feb 18, 2019 • edited Loading

simonbyrne commented Feb 18, 2019 • edited Loading

simonbyrne commented Feb 19, 2019

StefanKarpinski commented Feb 20, 2019

simonbyrne commented Feb 20, 2019

fp4code commented Mar 27, 2019 • edited Loading

cafaxo commented Feb 13, 2019 •

edited

Loading

stevengj commented Feb 14, 2019 •

edited

Loading

stevengj commented Feb 14, 2019 •

edited

Loading

stevengj commented Feb 14, 2019 •

edited

Loading

cafaxo commented Feb 14, 2019 •

edited

Loading

stevengj commented Feb 14, 2019 •

edited

Loading

cafaxo commented Feb 15, 2019 •

edited

Loading

cafaxo commented Feb 18, 2019 •

edited

Loading

simonbyrne commented Feb 18, 2019 •

edited

Loading

fp4code commented Mar 27, 2019 •

edited

Loading