conj(sqrt(complex(-5,0))) is not the same as sqrt(conj(complex(-5,0)))

[ViralBShah]

An example from page 15 of http://www.cs.berkeley.edu/~wkahan/JAVAhurt.pdf

julia> conj(sqrt(complex(-5,-0)))
0.0 - 2.23606797749979im

julia> sqrt(conj(complex(-5,-0)))
0.0 + 2.23606797749979im

This is also the case with log:

julia> conj(log(complex(-5,0)))
1.6094379124341003 - 3.141592653589793im

julia> log(conj(complex(-5,0)))
1.6094379124341003 + 3.141592653589793im

Octave gives the same result in both cases.

[jiahao]

This is not necessarily a bug, but rather the reflection of some choice of branch cut embedded in the complex sqrt function. It is easy to verify that the two different answers correspond to different square roots of -5.

julia>x= conj(sqrt(complex(-5,-0)))
0.0 - 2.23606797749979im

julia> conj(x^2)
-5.000000000000001 + 0.0im

julia> conj(x)^2
-5.000000000000001 + 0.0im

julia>x= sqrt(conj(complex(-5,-0)))
0.0 + 2.23606797749979im

julia> conj(x)^2
-5.000000000000001 - 0.0im

julia> conj(x^2)
-5.000000000000001 - 0.0im

[ViralBShah]

According to Kahan: http://people.freebsd.org/~das/kahan86branch.pdf

[ViralBShah]

Yes, they are different square roots of -5. I am not sure if the branch cut we have chosen is by design. I have changed the label to a decision rather than a bug.

[ViralBShah]

Should conj not change the signbit of 0 in the second case here?

julia> conj(complex(5,-0))
5 + 0im

julia> conj(complex(5,0))
5 + 0im

[jiahao]

Yes, that is clearly an error. The Kahan article blames this as the root cause, which would in effect determine the behavior at branch cuts.

[jiahao]

FWIW the first link you provide clearly favors the ISO C99 standard. Here is what it has to say on this topic (emph. mine):

7.3.3 Branch cuts
1 Some of the functions below hav e branch cuts, across which the function is
discontinuous. For implementations with a signed zero (including all IEC 60559
implementations) that follow the specifications of annex G, the sign of zero distinguishes
one side of a cut from another so the function is continuous (except for format
limitations) as the cut is approached from either side. For example, for the square root
function, which has a branch cut along the negative real axis, the top of the cut, with
imaginary part +0, maps to the positive imaginary axis, and the bottom of the cut, with
imaginary part −0, maps to the negative imaginary axis.
2 Implementations that do not support a signed zero (see annex F) cannot distinguish the
sides of branch cuts. These implementations shall map a cut so the function is continuous
as the cut is approached coming around the finite endpoint of the cut in a counter
clockwise direction. (Branch cuts for the functions specified here have just one finite
endpoint.) For example, for the square root function, coming counter clockwise around
the finite endpoint of the cut along the negative real axis approaches the cut from above,
so the cut maps to the positive imaginary axis

[jiahao]

julia> 0==-0
true

julia> bits(0)
"0000000000000000000000000000000000000000000000000000000000000000"

julia> bits(-0)
"0000000000000000000000000000000000000000000000000000000000000000"

julia> bits(-0.0e0)
"1000000000000000000000000000000000000000000000000000000000000000"

julia> bits(0.0e0)
"0000000000000000000000000000000000000000000000000000000000000000"

julia> 0.0e0 == -0.0e0
true

Erm.

[pao]

Literal 0 is an Int. Isn't everything else here floating-point?

[jiahao]

Oh, right, duh.

[jiahao]

Ah, so part of my confusion is that complex(-5,-0) is equivalent to complex(-5, 0) since they are both Complex{Int64}. This however doesn't quite address the original branch cut issue:

julia> conj(complex(5.0, -0.0))
5.0 + 0.0im

julia> conj(complex(5.0, +0.0))
5.0 - 0.0im

julia> conj(sqrt(complex(-5,+0)))
0.0 - 2.23606797749979im

julia> conj(sqrt(complex(-5,-0)))
0.0 - 2.23606797749979im

julia> conj(sqrt(complex(-5.0e0,-0.0e0)))
0.0 - 2.23606797749979im

julia> sqrt(conj(complex(-5.0e0,-0.0e0)))
0.0 + 2.23606797749979im

julia> sqrt(conj(complex(-5.0e0,+0.0e0)))
0.0 + 2.23606797749979im

julia> conj(sqrt(complex(-5.0e0,+0.0e0)))
0.0 - 2.23606797749979im

[andreasnoack]

There seem to two problems as @jiahao points out. There is a branch cut issue for the complex sqrt. log seems okay.

julia> conj(log(complex(-5.,-0.)))
1.6094379124341003 + 3.141592653589793im

julia> log(conj(complex(-5.,-0.)))
1.6094379124341003 + 3.141592653589793im

The second problems is how to deal with Complex{Int}. Maybe this is not so important because most often when you want to use sqrt or log the type is likely to be Complex{Float}.

[jiahao]

This behavior makes sense then from a complex analysis standpoint. The Riemann surface of complex log is well-behaved except at complex 0, where the branch cuts serve to chop up the Riemann surface into single-valued domains, whereas the square root has nontrivial singularities at the branch cut which require a decision to be made about which surface to follow at the branch cut.

[jiahao]

Going back and reading the ISO C99 standard, it seems clear to me that the behavior we see conforms to the branch cut decisions made in the standard. Complex{Float*} is in case 1 (with signed zeros) and Complex{Int*} is in case 2 (no signed zeros, use counterclockwise convention).

[simonbyrne]

I don't think it is matching up to the C99 standard for floats:

julia> sqrt(complex(-1.0,-0.0))
0.0 + 1.0im

julia> sqrt(complex(-1.0,-1e-100))
5.0e-101 - 1.0im

[JeffBezanson]

I think there is a comparison in there where there should be a copysign.

[JeffBezanson]

Fixed in d3c66e9.

[jiahao]

Currently in d3c66e9:

julia> sqrt(conj(complex(-5.0,-0.0)))
0.0 + 2.23606797749979im

julia> sqrt(conj(complex(-5.0,0.0)))
0.0 - 2.23606797749979im

julia> conj(sqrt(complex(-5,-0)))
0.0 - 2.23606797749979im

julia> conj(sqrt(complex(-5,0)))
0.0 - 2.23606797749979im

[StefanKarpinski]

Is that good or bad? Note that since there's isn't a -0 integer, the last two will always have the same behavior.

[jiahao]

The behavior @simonbyrne pointed out now appears to conform to the C99 standard.

Case 1: signed zeros, sqrt(x+0.0_im) is positive imaginary, sqrt(x-0.0_im) is negative imaginary.

julia> sqrt(complex(-1.0,+0.0))
0.0 + 1.0im

julia> sqrt(complex(-1.0,+1e-100))
5.0e-101 + 1.0im

julia> sqrt(complex(-1.0,-0.0))
0.0 - 1.0im

julia> sqrt(complex(-1.0,-1e-100))
5.0e-101 - 1.0im

Case 2: no signed zeros, sqrt(x+0*im) is positive imaginary.

julia> sqrt(complex(-1,0))
0.0 + 1.0im

Since sqrt(complex(-1,-0)) caused confusion above, this behavior may be a good thing to highlight in the documentation to clarify how Complex{Int} works.

[ViralBShah]

I missed all the action here after filing. With #2845, we should hopefully be able to visit all the branch cuts and be compliant with C99.

One less reason for William Kahan to be unhappy with us now. :-)