angle of complex zeros

[StefanKarpinski]

One of these things is not like the other ones:

julia> angle(0.0 + 0.0im)
0.0

julia> angle(0.0 - 0.0im)
-0.0

julia> angle(-0.0 + 0.0im)
0.0

julia> angle(-0.0 - 0.0im)
-3.141592653589793

Arguably, these should all be NaN.

[lobingera]

none of the above is incorrect (as we have 2pi periodicity) and for sure should not be NaN.

[StefanKarpinski]

The limit of angle(z) as z approaches zero can be any value from [-π, π] depending on where it's approaching from. That may be what you mean by saying that "none of the above incorrect" but they're also not correct either – we don't generally give random answers at complex singularities.

[StefanKarpinski]

I can see that defining angle(0) == 0 might be practical, but in that case all of the complex zeros should produce that result up to sign.

[lobingera]

I might be biased here, as i live happily with matlabs angle(complex(0,0)) =0 ... and some parts of my simulator survive singularities on that (antenna patterns with an NaN angle tend to be undefined).
It's a tricky thing, in doubt, i'd just follow the mainstream.

[GunnarFarneback]

One of these things is not like the other ones:

Yes,

julia> angle(-0.0 + 0.0im)
0.0

produces a result in the wrong quadrant.
(I'm aware that's not what you meant and angle(-0.0 - 0.0im) quite possibly hit the right quadrant by accident.)
This all falls back to atan2, which defers to libm.

[BenLauwens]

The argument of a complex number with modulus zero is undefined (see Wolfram).

[lobingera]

Wikipedia means: "The polar angle for the complex number 0 is indeterminate, but arbitrary choice of the angle 0 is common." https://en.wikipedia.org/wiki/Complex_number

0 is a valid arbitray choice.

[StefanKarpinski]

So there are two issues here:

• We want angle(±0.0 ± 0.0im) == ±0.0 for all complex zeros.
• We want the sign of the zero angle returned to equal the sign of the imaginary zero part.

And tests, of course. This is ready for action if anyone wants to tackle it.

[StefanKarpinski]

I'm marking this as a bug based on the following analysis:

• If anyone is relying on the magnitude of angle(±0.0 ± 0.0im) being zero then their code is wrong in the -0.0 - 0.0im case and this is a bug fix for them.
• If anyone is relying on the sign of the returned zero, then their code is wrong in the -0.0 + 0.0im case and this is a bug fix for them.

If they are not relying on either of these things, then their code will not be affected by a fix.

[stevengj]

The current behavior of angle is the same as the behavior of atan2. Isn't this standardized?

[stevengj]

Actually, I take that back, it doesn't match atan2 Nevermind, I forgot to reverse the argument order.

[stevengj]

Part of the confusion here is due to the fact that -0.0 + 0.0im doesn't preserve the sign of the zero in the real part:

julia> z = Complex(-0.0, 0.0)
-0.0 + 0.0im

julia> real(z)
-0.0

julia> angle(z)
3.141592653589793

julia> atan2(imag(z), real(z))
3.141592653589793

julia> w = -0.0 + 0.0im
0.0 + 0.0im

julia> real(w)
0.0

julia> angle(w)
0.0

This stems from the lack of an imaginary type (see also #5468 and #10000), since IEEE rules give -0.0 + 0.0 === 0.0.

[stevengj]

I removed the "bug" label and added the "decision" label. I think the correct decision is "leave as-is". This behavior of atan2 for signed zeros seems to be extremely widespread across many languages (since it follows the ISO C standard), and the behavior of angle follows (and is also standardized in ISO C as the carg function).

For example, in Python I also get:

>>> import cmath
>>> cmath.phase(complex(-0.0, 0.0))
3.141592653589793
>>> cmath.phase(complex(-0.0, -0.0))
-3.141592653589793
>>> cmath.phase(-0.0 + 0.0j)
0.0

Note that Python has the same problem with the lack of an imaginary type, so you need to call the complex constructor directly as in Julia to reliably get the desired sign of zero.

[BenLauwens]

In gfortran-6 compilation of the following code gives an error:

program angle
  implicit none
  print *, atan2(+0.0, +0.0)
  print *, atan2(+0.0, -0.0)
  print *, atan2(-0.0, -0.0)
  print *, atan2(-0.0, +0.0)
end

angle.f90:3:23:

   print *, atan2(+0.0, +0.0)
                       1
Error: If first argument of ATAN2 (1) is zero, then the second argument must not be zero
angle.f90:4:23:

   print *, atan2(+0.0, -0.0)
                       1
Error: If first argument of ATAN2 (1) is zero, then the second argument must not be zero
angle.f90:5:23:

   print *, atan2(-0.0, -0.0)
                       1
Error: If first argument of ATAN2 (1) is zero, then the second argument must not be zero
angle.f90:6:23:

   print *, atan2(-0.0, +0.0)
                       1
Error: If first argument of ATAN2 (1) is zero, then the second argument must not be zero

If both the real and the imaginary part are zero, the most logical thing we can do is to throw an indeterminate exception. It is up to the user to decide what the value of the argument should be depending on the application.

[BenLauwens]

The problem arises from the fact that we are evaluating a limit:
lim_{x\rightarrow 0, y\rightarrow} Arctan(y/x)
and it all depends on how x and y go to zero. The use of signed zeros is in fact the limit of x and y going to zero at the same rate:

• angle(Complex(+0,+0)) = \pi/4
• angle(Complex(+0,-0)) = -\pi/4
• angle(Complex(-0,+0)) = 3\pi/4
• angle(Complex(-0,-0)) = -3\pi/4
  For one non-signed zeros we have the following results:
• angle(Complex(0,+0)) = \pi/2
• angle(Complex(0,-0)) = -\pi/2
• angle(Complex(+0,0)) = 0
• angle(Complex(-0,0)) = \pi
  When both are non-signed zeros the limit does not exist.
  If we only want to convert a complex number from cartesian representation to exponential representation, it does not matter because the complex number is formed by the couple (x,y) or (r, \theta).
  I think the most important thing is to document the behaviour. Whatever choice we are making, someone will be unhappy ... so sticking to the definition of atan2 seems fine but if both x and y are zero the behaviour is undefined (see wikipedia)

[stevengj]

I'd much prefer to stay with the current behavior (and document it). When writing things in polar form, for example, it is annoying to have an exception thrown at 0+0im, even though any angle is valid there. Similarly I suspect that most programs would rather atan2 pick something at (0,0) rather than throw an exception, since any angle at (0,0) will usually be acceptable. Fortran's behavior for atan2 was probably set in the 60s or 70s, long before IEEE floating-point standards, and I'd much rather agree with the ISO C standard's choice of behavior on this.

[stevengj]

(I'm extremely suspicious of discarding the ISO C behavior of a basic mathematical function based on people's gut reactions in a github issue. If the ISO atan function returned NaN, that would be reasonable to translate to an exception in Julia, but not if it returns a numeric value.)

[GunnarFarneback]

I think this behaviour is what you really want to get out of angle for signed zeros.

julia> z2pol(z) = (abs(z), angle(z))
z2pol (generic function with 1 method)

julia> pol2z(p) = Complex(p[1] * cos(p[2]), p[1] * sin(p[2]))
pol2z (generic function with 1 method)

julia> pol2z(z2pol(Complex(0.0, 0.0)))
0.0 + 0.0im

julia> pol2z(z2pol(Complex(0.0, -0.0)))
0.0 - 0.0im

julia> pol2z(z2pol(Complex(-0.0, 0.0)))
-0.0 + 0.0im

julia> pol2z(z2pol(Complex(-0.0, -0.0)))
-0.0 - 0.0im

Unfortunately it doesn't hold for Float32 since π rounds the wrong way in that case.

julia> pol2z(z2pol(Complex(-0.0f0, 0.0f0)))
-0.0f0 - 0.0f0im

julia> pol2z(z2pol(Complex(-0.0f0, -0.0f0)))
-0.0f0 + 0.0f0im

[stevengj]

By the way, even thoughatan2(0,0) is an indeterminate form in the limit sense, it can still be convenient and appropriate to assign it a value in a floating-point library. For another example where we (and most other fp libraries) assign a value to an indeterminate form, see:

julia> 0.0^0.0
1.0

julia> (0.0 + 0.0im) ^ (0.0 + 0.0im)
1.0 + 0.0im

[stevengj]

See also the discussion of this when it came up in Octave; they decided to follow ISO/IEC. Matlab's atan2 returns 0 at the origin for all cases (with what sign I cannot tell...Matlab doesn't seem to give you access to the sign of zero), but this is probably a legacy issue: Matlab pre-dates IEEE 754 and ISO/IEC C, and it hides the sign of 0.0 from you too.

[stevengj]

The R atan2 function apparently follows ISO/IEC, though I couldn't find documentation of this. Ruby's atan2 function also follows ISO/IEC, though again I couldn't find documentation of this.

[simonbyrne]

Also, the behaviour of atan2 is also specified in the IEEE754 spec, and ISO 10967-2 (Language independent arithmetic), so I think deviating from this requires a pretty strong justification.

Another advantage of the current behaviour is that exp(log(complex(x,y))) === complex(x,y) for all signed zero x,y.

[simonbyrne]

On the other hand, -0.0 + 0.0im !== complex(-0.0,0.0) is perhaps a good argument for an imaginary type.

[simonbyrne]

@GunnarFarneback's comment is equivalent, but he is correct in that this doesn't hold for Float32.

julia> log(Complex(-0.0,0.0))
-Inf + 3.141592653589793im

julia> exp(log(Complex(-0.0,0.0)))
-0.0 + 0.0im

julia> log(Complex(-0f0,0f0))
-Inf32 + 3.1415927f0im

julia> exp(log(Complex(-0f0,0f0)))
-0.0f0 - 0.0f0im

which as he correctly points out, is due to

julia> Float32(pi) > pi
true

The IEEE754-2008 spec has the following to say on the matter:

For some formats under some rounding attributes the rounded magnitude range of atan (atan2) may exceed the unrounded magnitude of π/2 (π). In those cases, an anomalous manifold jump may occur under the inverse function for which the careful programmer should account.

However ISO 10967-2 says that arc_F (x, y) (which is their name for atan2, though they switch the order of the arguments) should give down_F (π) if x=−0 and y=0 (where down_F is the rounding function that rounds toward negative infinity).

So basically, the two specs disagree.

[simonbyrne]

For what it's worth, the OS X libm (unlike openlibm and glibc) appears to follow the ISO 10967 spec:

julia> ccall((:atan2f,:libm), Float32, (Float32, Float32), 0f0, -0f0) < pi
true

[simonbyrne]

And of course, this issue was discussed on the Octave lists as well

[simonbyrne]

I don't think we're realistically going to change the behaviour of atan2 (though if anyone finds out the reason for the discrepancy between the 2 specs, I would be grateful).

[StefanKarpinski]

Could we have a note in the docs for angle indicating what is done and why?

[ron-wolf]

@StefanKarpinski I’m down to write a pull request to add this to the docs. Is there any particular format notices like this are usually put in, or can it just be appended?

[simonbyrne]

Use !!! note block and indent, e.g. similar to

julia/base/floatfuncs.jl

Lines 96 to 116 in e58bc72

	!!! note
	Rounding to specified digits in bases other than 2 can be inexact when
	operating on binary floating point numbers. For example, the [`Float64`](@ref)
	value represented by `1.15` is actually *less* than 1.15, yet will be
	rounded to 1.2.
	# Examples
	```jldoctest; setup = :(using Printf)
	julia> x = 1.15
	1.15
	julia> @sprintf "%.20f" x
	"1.14999999999999991118"
	julia> x < 115//100
	true
	julia> round(x, digits=1)
	1.2
	```