Possible bug in unique/Set

I used unique on a vector of [nan, 0.0, -0.0, nan] as a surrogate, but also checked containers.Map in Matlab and Association in Mathematica, and std::map for c++. Here is a link to the comparisons I made with a little more detail. The link allows edits in case anyone wants to comment on that directly.

This might reveal a giant gap in my knowledge. I don’t understand the implications of user-extensible mechanisms and the difference between the two sets of comparative operators

That would go a long way. So for floating point numbers (not NaNs) the behaviors of <(a,b) and isless(a,b) would be consistent as would the behaviors of ==(a,b) and isequal(a,b). The functions and infix operators would still differ for NaNs, but that is OK by me.

I guess you are referring to breaking existing Julia code, since sorting NaNs differs so much across languages.

No need to change < operator wrt NaNs, They can remain non-comparable. So, if a list does not contain NaNs, then sorting and hashing using the infix operators, < and ==, would be totally ordered and consistent with the functional forms, isless and isequal. But, only the functional forms would have total ordering if NaNs were present, while the infix operators would have strict-weak ordering.

The only change that’s really plausible would be hashing and sorting -0.0 and 0.0 the same.

Unfortunately, I agree.

What about making min and max map to the processor intrinsics? Yes, the processor intrinsics vminsd and vmaxsd are special for NaN (they have a strong bias for returning the second argument).

Or at least make it branchfree? The following does afaik reproduce the current NaN behavior of min (also with the bias towards returning the second argument if both are NaN):

julia> fmin(x,y)=ifelse(x<y,x,y);
julia> gmin(x,y)=ifelse(x<y, x, ifelse(x>y, y, ifelse(y==y, x, y)));

julia> v=rand(1000).-0.5; w=rand(1000).-0.5; z=similar(v);
julia> using BenchmarkTools

julia> @btime broadcast!($min,$z,$v,$w);
  1.527 ��s (0 allocations: 0 bytes)
#06:   1.924 ��s (4 allocations: 96 bytes)

julia> @btime broadcast!($fmin,$z,$v,$w);
  210.393 ns (0 allocations: 0 bytes)
#06:   448.500 ns (4 allocations: 96 bytes)
julia> @btime broadcast!($gmin,$z,$v,$w);
  689.187 ns (0 allocations: 0 bytes)
#06:   873.531 ns (4 allocations: 96 bytes)

julia> v=rand(Float32,1000).-Float32(0.5); w=rand(Float32,1000).-Float32(0.5); z=similar(v);

julia> @btime broadcast!($min,$z,$v,$w);
  794.474 ns (0 allocations: 0 bytes)
#0.6  1.118 ��s (4 allocations: 96 bytes)
julia> @btime broadcast!($fmin,$z,$v,$w);
  129.841 ns (0 allocations: 0 bytes)
 #0.6 370.486 ns (4 allocations: 96 bytes)
julia> @btime broadcast!($gmin,$z,$v,$w);
  364.726 ns (0 allocations: 0 bytes)
#0.6  604.028 ns (4 allocations: 96 bytes)

Perhaps a better approach would be having @fast_math versions of them that do whatever the processor does. I do wish the instructions sets put a little more thought into what would be a pleasant high-level behavior.

1 Like

Because they don’t make any sense! See https://github.com/JuliaLang/julia/issues/7866 for the full discussion.

Issue opened: https://github.com/JuliaLang/julia/issues/27545.

If a data file contains a -0 that needs to be interpreted as a floating point negative zero, that has to be handled by the code parsing the data file. -0 as julia syntax is a separate issue; it’s an expression that evaluates to 0, the same way 1+1 evaluates to 2.

Very interesting comparison. Thanks for doing that and sharing it!

In systems where you either can’t define new user types or you can’t extend existing operators to new types, you don’t need a name for the relations used for hashing and sorting. In your analysis, C++, R and Julia all have different total orders for some things than others—the normal </== order and what they use for sorting and hashing. Julia gives these different orderings names (==/< and isequal/isless) and lets you extend them. The fact that they’re implicit in other languages and not extensible or accessible to users doesn’t mean they don’t exist.

Mathematica and Matlab, on the other hand are consistent and just use < and == everywhere, so they truly only have the one ordering. But they both have to sacrifice something to accomplish this. Mathematica sacrifices IEEE 754 conformance. At this point we’re not going to do that.

What Matlab sacrifices is that it considers NaN distinct from itself even in uniquing and hashing. That may not seem like a big deal, but it means that if some long running process ends up using NaN as a dictionary key repeatedly, it will use an arbitrary amount of memory. For example, this makes it trivially easy to create a denial of service attack if you can convince a program to use a floating-point value that you can inject: just generate NaN over and over. Of course, that’s not a huge concern because only a crazy person would do something like deploy Matlab on the Internet (confession: I am such a crazy person; it was a long time ago), but in a language that you can reasonably do something like building a network service in—which includes Julia—that’s really not ok. There’s also a somewhat more philosophical argument to be made that every equivalence relation should be reflexive and that they should never be finer grained that egal. Bottom line: the same bit pattern should definitely be equal to itself and hash the same as itself.

Sorting is not really the issue—no one really cares about the sort order of NaNs. The issue is changing x ≤ NaN from always false to always true. That means any code that checks, for example, 0.0 < x and then does something based on that check, assuming that the value x will be a positive finite value or Inf, will suddenly be broken because x could also be NaN. At this point, that’s a lot of existing Julia code and a lot of code that’s waiting to be ported from Fortran, C, R, Matlab, Python or any other language that conforms to IEEE 754 with respect to NaNs. Basically all of them besides Mathematica.

7 Likes

About NaN:

There are two flavors of NaN, signaling and quiet. Signaling NaNs are “abort immediately” indicators and Julia does not deal with them at present. LLVM has a predicate isSignaling that is the partner of isNaN which would be used to distinguish them. isNaN returns true for either flavor of NaN. Quiet NaNs are “propogate self” tokens that will percolate from wherever they are introduced through to the end of the encompassing calculation chain. Put another way, op(x, NaN) === ope(NaN, x) === NaN.

There are many distinct quiet NaNs. Julia (and most other languages) use the NaN with a payload (an integer that is carried within the NaN) of zero as what one gets writing NaN. We can set the payload as we like, and the NaN with that payload will propagate. This is useful when one wants to propogate information along with the occurance of an inapproprate or illegal value.

julia> reinterpret(UInt64, NaN), reinterpret(UInt32, NaN32), reinterpret(UInt16, NaN16)
(0x7ff8000000000000, 0x7fc00000, 0x7e00)


julia> using QNaNs

julia> qnan(Int64(1)), qnan(Int32(-1)), qnan(Int16(3))
(NaN, NaN32, NaN16)

julia> reinterpret(UInt64, ans[1]), reinterpret(UInt32, ans[2]), reinterpret(UInt16, ans[3])
(0x7ff8000000000001, 0xffc00001, 0x7e03)

julia> qnan(Int64(1)), qnan(Int32(-1)), qnan(Int16(3))
(NaN, NaN32, NaN16)

julia> qnan.(ans)
(1, -1, 3)
1 Like

“Broken” may be too harsh. The code would work, in that the NaN would propogate, as it should. What would be lost is any special handling of NaNs that might appear in the else part of the code. In most code that does check for NaN, that check is first, either explicitly isnan or implicitly if !isfinite.

Also, consider what happens when a function checks for NaN. Typically, it throws an exception – or lets the NaN go forward. None of these points contradicts yours; it may soften though.

The Rust solution is also interesting: they deal with NaN by disallowing certain comparisons on floats.

https://www.rust-lang.org/en-US/faq.html#why-cant-i-compare-floats

I find it a really nice solution, actually. “Hashing on Float64? Just say no!” :wink:

1 Like