NaN vs wild (or, what's a DomainError, really?)

[jiahao]

(Context: this issue has come up recently in #4967 and elsewhere)

We currently don't treat the results of indeterminate computations consistently between real and complex arithmetic.

Example 1: cosine

julia> cos(Inf)
ERROR: DomainError
 in cos at math.jl:277

julia> cos(complex(Inf,0))
NaN - 0.0im

Example 2: inverse

inverse of 0

julia> inv(0.0)
Inf

julia> inv(complex(0.0, 0.0))
NaN + NaN*im

This illustrates a problem with infinity and signed zero. Unlike real infinities, of which there are only two FloatingPoint values, you can have 13 possible representations of complex infinities which convey different information about the phase (argument) of the complex infinity:

complex(+Inf, +0.0) #phase is exactly 0 or approaches 0 from the first quadrant
complex(+Inf, +Inf) #phase is in (0, pi/2), i.e. the first quadrant excluding the edges
complex(+0.0, +Inf) #phase is exactly pi/2 or approaches pi/2 from the first quadrant
complex(-0.0, +Inf) #phase approaches pi/2 from the second quadrant
complex(-Inf, +Inf) #phase is in (pi/2, pi), i.e. the second quadrant excluding the edges
complex(-Inf, +0.0) #phase is exactly pi or approaches pi from the second quadrant
complex(-Inf, -0.0) #phase approaches pi from the third quadrant
complex(-Inf, -Inf) #phase is in (pi, 3pi/2), i.e. the third quadrant excluding the edges
complex(-0.0, -Inf) #phase approaches 3pi/2 from the third quadrant
complex(+0.0, -Inf) #phase is exactly 3pi/2 or approaches 3pi/2 from the fourth quadrant
complex(+Inf, -Inf) #phase is in (3pi/2, 2pi), i.e. the fourth quadrant excluding the edges
complex(+Inf, -0.0) #phase approaches 2pi from the fourth quadrant
complex(NaN, NaN) #phase cannot be determined to lie in exactly one of the above regions
#                  and hence the infinity has no valid non-NaN representation in floating point

The result is correct if we work with unsigned zeros. However, after accounting for the signed zeros in complex(+0.0, +0.0), this result should be complex(+Inf, +Inf). A DivideError might also be a reasonable alternative here.

inverse of infinities

The inverse mapping is also problematic:

julia> inv(Inf)
0.0

julia> inv(complex(Inf,0))
0.0 - 0.0im

julia> inv(complex(Inf,Inf))
NaN + NaN*im

inverse of indeterminates

julia> inv(NaN)
NaN

julia> inv(complex(0,NaN))
NaN + NaN*im

julia> inv(complex(NaN, 0))
NaN + NaN*im

julia> inv(complex(NaN, NaN))
NaN + NaN*im

Is it meaningful to distinguish between these three possible complex NaNs? It seems silly in this example, but consider also:

julia> complex(0, NaN) + complex(0, NaN)
complex(0.0,NaN)

julia> complex(0, NaN) * complex(0, NaN)
NaN + NaN*im

etc.

Example 3: roots (nonintegral powers)

julia> sqrt(-1)
ERROR: DomainError

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

julia> (-1)^(-1/2)
NaN

julia> complex(-1)^(-1/2)
6.123233995736766e-17 - 1.0im

Whereas sqrt(-1) returns the notorious DomainError, the inverse square root (as computed by x->x^-1/2 does not, but returns a NaN instead. This use of NaN is sanctioned by IEEE 754 in the specific case of a real operation with no real output.

tl;dr: if a floating-point computation returns an indeterminate value, when should it return a NaN (or any of its complex variants), and when should it throw an error like DivideError or DomainError? By my reading of IEEE 754, both of these behaviors are allowed (throwing an error would correspond to a signaling NaN which is trapped by the error handler). Each of these examples in isolation show valid behavior; however, we should be consistent.

[StefanKarpinski]

I think we should throw an error whenever the condition was explicitly checked for anyway.

[jiahao]

To summarize the prevailing thinking with @StefanKarpinski @loladiro @JeffBezanson:

Operations whose results would produce NaNs should throw Exceptions if the cost of testing a NaN would not dominate the cost of the operation.

By this rule of thumb, elementary arithmetic on real floating-point numbers would be allowed to return NaNs if necessary, but more complicated operations like complex division should throw a DivideError.

[ivarne]

I don't understand. What operations should be allowed to return NaN, (unless the argument is NaN)? Only /(a::FloatingPoint, b::FloatingPoint)?

[jiahao]

I think we only have a clear idea that real elementary arithmetic (+, -, x, /) should be allowed to produce and propagate NaNs, since the cost of testing for a NaN or inputs that would produce NaNs would be comparable (or more expensive) than the actual operation itself, but not much else.

[ivarne]

So all other functions (also those who does not return NaN for any real input) must check for NaN on input and throw DomainError?

[JeffBezanson]

No, NaN input can just propagate. We mostly don't want to produce NaNs from non-NaNs.

[StefanKarpinski]

The poster child is complex division: rather than generating a complex NaN when dividing by zero, just raise an error. Since complex division is so, um, complex, checking for a zero divisor is not exactly going to affect it's performance in any meaningful way. Raising and error is simpler and alerts the user to the problem sooner.

[nalimilan]

Why not, but if that's the only operation which would be affected, is it worth introducing an inconsistency? It would be good to have a simple rule for which functions return NaN, and which throw an error (I mean, an user-understandable rule, not something based on computing costs).

[JeffBezanson]

My guess at the moment is that the rule is: only real floating-point +, -, *, / can produce NaNs from non-NaNs.

[JeffBezanson]

I would also add complex floating-point +, -, *, / to that list. The 1/complex(0.0,0.0) case is a divide-by-zero error, so you can still get NaNs from complex / in other ways.