implement extrema with mapreduce machinery

Update multidimensional.jl

Author: N5N3

base/multidimensional.jl

@@ -1746,41 +1746,52 @@ of `A`.
     This method requires Julia 1.2 or later.
 """
 extrema(f, A::AbstractArray; dims=:) = _extrema_dims(f, A, dims)
-
-_extrema_dims(f, A::AbstractArray, ::Colon) = _extrema_itr(f, A)
-
-function _extrema_dims(f, A::AbstractArray, dims)
-    sz = size(A)
-    for d in dims
-        sz = setindex(sz, 1, d)
-    end
-    T = promote_op(f, eltype(A))
-    B = Array{Tuple{T,T}}(undef, sz...)
-    return extrema!(f, B, A)
-end
-
-@noinline function extrema!(f, B, A)
-    require_one_based_indexing(B, A)
-    sA = size(A)
-    sB = size(B)
-    for I in CartesianIndices(sB)
-        fAI = f(A[I])
-        B[I] = (fAI, fAI)
-    end
-    Bmax = CartesianIndex(sB)
-    @inbounds @simd for I in CartesianIndices(sA)
-        J = min(Bmax,I)
-        BJ = B[J]
-        fAI = f(A[I])
-        if fAI < BJ[1]
-            B[J] = (fAI, BJ[2])
-        elseif fAI > BJ[2]
-            B[J] = (BJ[1], fAI)
-        end
+_extrema_dims(f, A::AbstractArray, dims) = mapreduce(x -> (fx = f(x); (fx, fx)), _extrema_op, A; dims)
+extrema!(B, A) = extrema!(identity, B, A)
+extrema!(f, B, A) = mapreduce!(x -> (fx = f(x); (fx, fx)), _extrema_op, B, A)
+_extrema_op((a, b), (c, d)) = min(a, c), max(b, d)
+# TODO: fuse this with min/max
+_isless(x::T, y::T) where {T<:Union{Float32,Float64}}= signbit(x - y)
+_isless(x::T, y::T) where {T<:AbstractFloat}= (x < y) | (signbit(x) > signbit(y))
+function _extrema_op(x::NTuple{2,T}, y::NTuple{2,T}) where {T<:AbstractFloat}
+    (x1, x2), (y1, y2) = x, y
+    z1 = ifelse(isnan(x1) | ~isnan(y1) & _isless(x1, y1), x1, y1)
+    z2 = ifelse(isnan(x1) | ~isnan(y1) & _isless(y2, x2), x2, y2)
+    z1, z2
+end
+function mapreduce_impl(f, ::typeof(_extrema_op), A::AbstractArrayOrBroadcasted, fi::Int, la::Int)
+    @inline elf(i) = @inbounds f(A[i])
+    Eltype = _return_type(elf, Tuple{Int})
+    Eltype <: NTuple{2,AbstractFloat} ||
+        return invoke(mapreduce_impl,Tuple{Any,Any,AbstractArrayOrBroadcasted,Int,Int},f,_extrema_op,A,fi,la)
+    v1 = v2 = v3 = v4 = Eltype((Inf,-Inf))
+    iter = fi:4:la-4
+    for i in iter
+        v1 = _extrema_op(v1, elf(i))
+        v2 = _extrema_op(v2, elf(i+1))
+        v3 = _extrema_op(v3, elf(i+2))
+        v4 = _extrema_op(v4, elf(i+3))
+    end
+    v = _extrema_op(_extrema_op(v1, v2), _extrema_op(v3, v4))
+    for i in last(iter)+4:la
+        v = _extrema_op(v, elf(i))
+    end
+    v
+end
+function reducedim_init(f, ::typeof(_extrema_op), A::AbstractArray, region)
+    ri = reduced_indices(A, region)
+    any(i -> isempty(axes(A, i)), region) && _empty_reduce_error()
+    A1 = view(A, ri...)
+    IT = eltype(A)
+    if missing isa IT
+        RT = promote_typejoin_union(_return_type(i -> f(i)[1], Tuple{nonmissingtype(IT)}))
+        T = Union{Tuple{RT,RT},Tuple{Missing,Missing}}
+    else
+        RT = promote_typejoin_union(_return_type(i -> f(i)[1], Tuple{IT}))
+        T = Union{Tuple{RT,RT}}
     end
-    return B
+    map!(f, reducedim_initarray(A,region,undef,T), A1)
 end
-extrema!(B, A) = extrema!(identity, B, A)
 
 # Show for pairs() with Cartesian indices. Needs to be here rather than show.jl for bootstrap order
 function Base.showarg(io::IO, r::Iterators.Pairs{<:Integer, <:Any, <:Any, T}, toplevel) where T <: Union{AbstractVector, Tuple}