optimize min/max and related for floats

[StefanKarpinski]

Someone on Discourse recently posted about findmax on floats being 4x slower than the C version. We do a bunch of special handling that doesn't necessarily match what the fast native instructions do, but I did some poking around and I think we could do much better than we currently are:

• On ARM the fmax and fmin instructions already do the right thing (propagate NaNs), so we should make sure that we emit the simple code there.
• On x64, it seems like the native instructions actually happen to implement isless (all NaNs sort after all non-NaNs), and we should make sure that this is used where appropriate.
• For maximum and findmax this is actually already the right order! So we should make sure that we just use the native instruction there.
• For minimum and findmin this is not the right order, but we can just negate, take the max, and then negate again at the end and get the right answer.

Am I missing something here, or could we be doing much better?

[giordano]

I think #56371 already moved things in the right direction?

[oscardssmith]

looks like 1.12 is ~6% faster here... The LLVM and assembly still seem fairly bad though.

julia> a = Float64.(1:128);
julia> @btime findmax($a)
  127.487 ns (0 allocations: 0 bytes) # 1.11
  120.244 ns (0 allocations: 0 bytes) # 1.12

julia> function naive_findmax(w)
           x = @fastmath foldl(max, w)
           i = findfirst(==(x), w)::Int
           x, i
       end
naive_findmax (generic function with 6 methods)

julia> @btime naive_findmax($a)
  72.061 ns (0 allocations: 0 bytes) # 1.11/1.12
(128.0, 128)

[mikmoore]

This issue overlaps with the ambitions of #45581, although that PR is abandoned because of the headaches of dimensional reduction, my personal anxiety over the implementation complexity, and the hope that maybe someday(tm) we'd get nice x86 intrinsics (maybe the recent LLVM intrinsics already help a tad?).

Using that PR as a template yields this one-pass function with NaN and signed-zero handling:

function _fpord(x::T) where T<:Base.IEEEFloat
    # converts IEEEFloat to BitSigned with ordering [-NaN, -Inf, ..., -0.0, +0.0, ..., +Inf, +NaN]
    # note: the precise implementation of this function can yield significantly different performance
    #    this version may not be the best
    u = reinterpret(Base.inttype(T), x)
    uf = xor(u, typemax(u)) # u with all non-sign bits flipped
    o = ifelse(signbit(u), uf, u)
    return o
end

function fpord_max(x::T) where T<:Base.IEEEFloat
    o = _fpord(x) # ordering [-NaN, -Inf, ..., -0.0, +0.0, ..., +Inf, +NaN]
    omin = _fpord(typemin(T)) # smallest representable value of x
    return o + (typemin(omin) - omin) # shift so that typemin is smallest
    # ordering [-Inf, ... -0.0, +0.0, ... +Inf, +NaN, -NaN]
end

function findmax2(x)
    o,i = mapreduce(ix -> (fpord_max(last(ix)), first(ix)), (a,b) -> first(b)>first(a) ? b : a, pairs(x))
    return x[i],i # or could reverse-transform o
end

And now a comparison (v1.11):

julia> a = Float64.(1:128);

julia> @btime findmax($a)
  112.501 ns (0 allocations: 0 bytes)
(128.0, 128)

julia> @btime naive_findmax($a)
  68.532 ns (0 allocations: 0 bytes)
(128.0, 128)

julia> @btime findmax2($a)
  71.915 ns (0 allocations: 0 bytes)
(128.0, 128)

A longer input does favor naive_findmax, though recall that it has the benefit of mishandling zeros and NaNs:

julia> b = randn(10^5); b[end-10] = 1e3;

julia> @btime findmax($b)
  82.235 μs (0 allocations: 0 bytes)
(1000.0, 99990)

julia> @btime naive_findmax($b)
  38.068 μs (0 allocations: 0 bytes)
(1000.0, 99990)

julia> @btime findmax2($b)
  50.621 μs (0 allocations: 0 bytes)
(1000.0, 99990)

[mikmoore]

On x64, it seems like the native instructions actually happen to implement isless (all NaNs sort after all non-NaNs), and we should make sure that this is used where appropriate.

This isn't true. As pointed out on Discourse (after this issue was open, in fairness), the x86 intrinsics are equivalent to the common idioms x < y ? x : y for min and x > y ? x : y for max. This means that they mishandle (for the purposes of Julia semantics) NaNs and signed zeros.