diff --git a/NEWS.md b/NEWS.md
index 3f23aaf573b72..98baf83ad1c1e 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -99,6 +99,7 @@ Standard library changes
   arithmetic to error if the result may be wrapping. Or use a package such as SaferIntegers.jl when
   constructing the range. ([#40382])
 * TCP socket objects now expose `closewrite` functionality and support half-open mode usage ([#40783]).
+* `extrema` now supports `init` keyword argument ([#36265], [#43604]).
 * Intersect returns a result with the eltype of the type-promoted eltypes of the two inputs ([#41769]).
 * `Iterators.countfrom` now accepts any type that defines `+`. ([#37747])
 
diff --git a/base/compiler/compiler.jl b/base/compiler/compiler.jl
index c265512afcbf6..aafd80c1127ff 100644
--- a/base/compiler/compiler.jl
+++ b/base/compiler/compiler.jl
@@ -131,6 +131,19 @@ include("compiler/abstractinterpretation.jl")
 include("compiler/typeinfer.jl")
 include("compiler/optimize.jl") # TODO: break this up further + extract utilities
 
+# required for bootstrap
+# TODO: find why this is needed and remove it.
+function extrema(x::Array)
+    isempty(x) && throw(ArgumentError("collection must be non-empty"))
+    vmin = vmax = x[1]
+    for i in 2:length(x)
+        xi = x[i]
+        vmax = max(vmax, xi)
+        vmin = min(vmin, xi)
+    end
+    return vmin, vmax
+end
+
 include("compiler/bootstrap.jl")
 ccall(:jl_set_typeinf_func, Cvoid, (Any,), typeinf_ext_toplevel)
 
diff --git a/base/exports.jl b/base/exports.jl
index 4daddd664a21a..4d6ceef6a13c0 100644
--- a/base/exports.jl
+++ b/base/exports.jl
@@ -385,6 +385,7 @@ export
     eachindex,
     eachrow,
     eachslice,
+    extrema!,
     extrema,
     fill!,
     fill,
diff --git a/base/multidimensional.jl b/base/multidimensional.jl
index 0843cec35c26e..74ef782d92c69 100644
--- a/base/multidimensional.jl
+++ b/base/multidimensional.jl
@@ -1695,80 +1695,6 @@ _unique_dims(A::AbstractArray, dims::Colon) = invoke(unique, Tuple{Any}, A)
     end
 end
 
-"""
-    extrema(A::AbstractArray; dims) -> Array{Tuple}
-
-Compute the minimum and maximum elements of an array over the given dimensions.
-
-# Examples
-```jldoctest
-julia> A = reshape(Vector(1:2:16), (2,2,2))
-2×2×2 Array{Int64, 3}:
-[:, :, 1] =
- 1  5
- 3  7
-
-[:, :, 2] =
-  9  13
- 11  15
-
-julia> extrema(A, dims = (1,2))
-1×1×2 Array{Tuple{Int64, Int64}, 3}:
-[:, :, 1] =
- (1, 7)
-
-[:, :, 2] =
- (9, 15)
-```
-"""
-extrema(A::AbstractArray; dims = :) = _extrema_dims(identity, A, dims)
-
-"""
-    extrema(f, A::AbstractArray; dims) -> Array{Tuple}
-
-Compute the minimum and maximum of `f` applied to each element in the given dimensions
-of `A`.
-
-!!! compat "Julia 1.2"
-    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
-    end
-    return B
-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}
     print(io, "pairs(::$T)")
diff --git a/base/operators.jl b/base/operators.jl
index 17b5c739d6b38..52c64bb56c61e 100644
--- a/base/operators.jl
+++ b/base/operators.jl
@@ -504,58 +504,6 @@ julia> minmax('c','b')
 """
 minmax(x,y) = isless(y, x) ? (y, x) : (x, y)
 
-"""
-    extrema(itr) -> Tuple
-
-Compute both the minimum and maximum element in a single pass, and return them as a 2-tuple.
-
-# Examples
-```jldoctest
-julia> extrema(2:10)
-(2, 10)
-
-julia> extrema([9,pi,4.5])
-(3.141592653589793, 9.0)
-```
-"""
-extrema(itr) = _extrema_itr(identity, itr)
-
-"""
-    extrema(f, itr) -> Tuple
-
-Compute both the minimum and maximum of `f` applied to each element in `itr` and return
-them as a 2-tuple. Only one pass is made over `itr`.
-
-!!! compat "Julia 1.2"
-    This method requires Julia 1.2 or later.
-
-# Examples
-```jldoctest
-julia> extrema(sin, 0:π)
-(0.0, 0.9092974268256817)
-```
-"""
-extrema(f, itr) = _extrema_itr(f, itr)
-
-function _extrema_itr(f, itr)
-    y = iterate(itr)
-    y === nothing && throw(ArgumentError("collection must be non-empty"))
-    (v, s) = y
-    vmin = vmax = f(v)
-    while true
-        y = iterate(itr, s)
-        y === nothing && break
-        (x, s) = y
-        fx = f(x)
-        vmax = max(fx, vmax)
-        vmin = min(fx, vmin)
-    end
-    return (vmin, vmax)
-end
-
-extrema(x::Real) = (x, x)
-extrema(f, x::Real) = (y = f(x); (y, y))
-
 ## definitions providing basic traits of arithmetic operators ##
 
 """
diff --git a/base/reduce.jl b/base/reduce.jl
index 0581a72fb2862..13e1b69c79ede 100644
--- a/base/reduce.jl
+++ b/base/reduce.jl
@@ -604,7 +604,7 @@ julia> prod(1:5; init = 1.0)
 """
 prod(a; kw...) = mapreduce(identity, mul_prod, a; kw...)
 
-## maximum & minimum
+## maximum, minimum, & extrema
 _fast(::typeof(min),x,y) = min(x,y)
 _fast(::typeof(max),x,y) = max(x,y)
 function _fast(::typeof(max), x::AbstractFloat, y::AbstractFloat)
@@ -634,11 +634,6 @@ function mapreduce_impl(f, op::Union{typeof(max), typeof(min)},
     start = first + 1
     simdstop  = start + chunk_len - 4
     while simdstop <= last - 3
-        # short circuit in case of NaN or missing
-        (v1 == v1) === true || return v1
-        (v2 == v2) === true || return v2
-        (v3 == v3) === true || return v3
-        (v4 == v4) === true || return v4
         @inbounds for i in start:4:simdstop
             v1 = _fast(op, v1, f(A[i+0]))
             v2 = _fast(op, v2, f(A[i+1]))
@@ -785,6 +780,76 @@ Inf
 """
 minimum(a; kw...) = mapreduce(identity, min, a; kw...)
 
+"""
+    extrema(itr; [init]) -> (mn, mx)
+
+Compute both the minimum `mn` and maximum `mx` element in a single pass, and return them
+as a 2-tuple.
+
+The value returned for empty `itr` can be specified by `init`. It must be a 2-tuple whose
+first and second elements are neutral elements for `min` and `max` respectively
+(i.e. which are greater/less than or equal to any other element). As a consequence, when
+`itr` is empty the returned `(mn, mx)` tuple will satisfy `mn ≥ mx`. When `init` is
+specified it may be used even for non-empty `itr`.
+
+!!! compat "Julia 1.8"
+    Keyword argument `init` requires Julia 1.8 or later.
+
+# Examples
+```jldoctest
+julia> extrema(2:10)
+(2, 10)
+
+julia> extrema([9,pi,4.5])
+(3.141592653589793, 9.0)
+
+julia> extrema([]; init = (Inf, -Inf))
+(Inf, -Inf)
+```
+"""
+extrema(itr; kw...) = extrema(identity, itr; kw...)
+
+"""
+    extrema(f, itr; [init]) -> (mn, mx)
+
+Compute both the minimum `mn` and maximum `mx` of `f` applied to each element in `itr` and
+return them as a 2-tuple. Only one pass is made over `itr`.
+
+The value returned for empty `itr` can be specified by `init`. It must be a 2-tuple whose
+first and second elements are neutral elements for `min` and `max` respectively
+(i.e. which are greater/less than or equal to any other element). It is used for non-empty
+collections. Note: it implies that, for empty `itr`, the returned value `(mn, mx)` satisfies
+`mn ≥ mx` even though for non-empty `itr` it  satisfies `mn ≤ mx`.  This is a "paradoxical"
+but yet expected result.
+
+!!! compat "Julia 1.2"
+    This method requires Julia 1.2 or later.
+
+!!! compat "Julia 1.8"
+    Keyword argument `init` requires Julia 1.8 or later.
+
+# Examples
+```jldoctest
+julia> extrema(sin, 0:π)
+(0.0, 0.9092974268256817)
+
+julia> extrema(sin, Real[]; init = (1.0, -1.0))  # good, since -1 ≤ sin(::Real) ≤ 1
+(1.0, -1.0)
+```
+"""
+extrema(f, itr; kw...) = mapreduce(ExtremaMap(f), _extrema_rf, itr; kw...)
+
+# Not using closure since `extrema(type, itr)` is a very likely use-case and it's better
+# to avoid type-instability (#23618).
+struct ExtremaMap{F} <: Function
+    f::F
+end
+ExtremaMap(::Type{T}) where {T} = ExtremaMap{Type{T}}(T)
+@inline (f::ExtremaMap)(x) = (y = f.f(x); (y, y))
+
+# TODO: optimize for inputs <: AbstractFloat
+@inline _extrema_rf((min1, max1), (min2, max2)) = (min(min1, min2), max(max1, max2))
+
 ## findmax, findmin, argmax & argmin
 
 """
diff --git a/base/reducedim.jl b/base/reducedim.jl
index 02b2345038bcc..d55db2768e62b 100644
--- a/base/reducedim.jl
+++ b/base/reducedim.jl
@@ -77,15 +77,14 @@ end
 ## initialization
 # initarray! is only called by sum!, prod!, etc.
 for (Op, initfun) in ((:(typeof(add_sum)), :zero), (:(typeof(mul_prod)), :one))
-    @eval initarray!(a::AbstractArray{T}, ::$(Op), init::Bool, src::AbstractArray) where {T} = (init && fill!(a, $(initfun)(T)); a)
+    @eval initarray!(a::AbstractArray{T}, ::Any, ::$(Op), init::Bool, src::AbstractArray) where {T} = (init && fill!(a, $(initfun)(T)); a)
 end
 
-for Op in (:(typeof(max)), :(typeof(min)))
-    @eval initarray!(a::AbstractArray{T}, ::$(Op), init::Bool, src::AbstractArray) where {T} = (init && copyfirst!(a, src); a)
-end
+initarray!(a::AbstractArray{T}, f, ::Union{typeof(min),typeof(max),typeof(_extrema_rf)},
+    init::Bool, src::AbstractArray) where {T} = (init && mapfirst!(f, a, src); a)
 
 for (Op, initval) in ((:(typeof(&)), true), (:(typeof(|)), false))
-    @eval initarray!(a::AbstractArray, ::$(Op), init::Bool, src::AbstractArray) = (init && fill!(a, $initval); a)
+    @eval initarray!(a::AbstractArray, ::Any, ::$(Op), init::Bool, src::AbstractArray) = (init && fill!(a, $initval); a)
 end
 
 # reducedim_initarray is called by
@@ -165,6 +164,42 @@ for (f1, f2, initval, typeextreme) in ((:min, :max, :Inf, :typemax), (:max, :min
         end
     end
 end
+
+function reducedim_init(f::ExtremaMap, op::typeof(_extrema_rf), A::AbstractArray, region)
+    # First compute the reduce indices. This will throw an ArgumentError
+    # if any region is invalid
+    ri = reduced_indices(A, region)
+
+    # Next, throw if reduction is over a region with length zero
+    any(i -> isempty(axes(A, i)), region) && _empty_reduce_error()
+
+    # Make a view of the first slice of the region
+    A1 = view(A, ri...)
+
+    isempty(A1) && return map(f, A1)
+    # use the max/min of the first slice as initial value for non-empty cases
+    v0 = reverse(mapreduce(f, op, A1)) # turn minmax to maxmin
+
+    T = _realtype(f.f, promote_union(eltype(A)))
+    Tmin = v0[1] isa T ? T : typeof(v0[1])
+    Tmax = v0[2] isa T ? T : typeof(v0[2])
+
+    # but NaNs and missing need to be avoided as initial values
+    if v0[1] isa Number && isnan(v0[1])
+        v0 = oftype(v0[1], Inf), oftype(v0[2], -Inf)
+    elseif isunordered(v0[1])
+        # v0 is missing or a third-party unordered value
+        # TODO: Some types, like BigInt, don't support typemin/typemax.
+        # So a Matrix{Union{BigInt, Missing}} can still error here.
+        v0 = typemax(nonmissingtype(Tmin)), typemin(nonmissingtype(Tmax))
+    end
+    # v0 may have changed type.
+    Tmin = v0[1] isa T ? T : typeof(v0[1])
+    Tmax = v0[2] isa T ? T : typeof(v0[2])
+
+    return reducedim_initarray(A, region, v0, Tuple{Tmin,Tmax})
+end
+
 reducedim_init(f::Union{typeof(abs),typeof(abs2)}, op::typeof(max), A::AbstractArray{T}, region) where {T} =
     reducedim_initarray(A, region, zero(f(zero(T))), _realtype(f, T))
 
@@ -435,7 +470,7 @@ julia> count!(<=(2), [1; 1], A)
 """
 count!(r::AbstractArray, A::AbstractArrayOrBroadcasted; init::Bool=true) = count!(identity, r, A; init=init)
 count!(f, r::AbstractArray, A::AbstractArrayOrBroadcasted; init::Bool=true) =
-    mapreducedim!(_bool(f), add_sum, initarray!(r, add_sum, init, A), A)
+    mapreducedim!(_bool(f), add_sum, initarray!(r, f, add_sum, init, A), A)
 
 """
     sum(A::AbstractArray; dims)
@@ -737,6 +772,74 @@ julia> minimum!([1 1], A)
 """
 minimum!(r, A)
 
+"""
+    extrema(A::AbstractArray; dims) -> Array{Tuple}
+
+Compute the minimum and maximum elements of an array over the given dimensions.
+
+See also: [`minimum`](@ref), [`maximum`](@ref), [`extrema!`](@ref).
+
+# Examples
+```jldoctest
+julia> A = reshape(Vector(1:2:16), (2,2,2))
+2×2×2 Array{Int64, 3}:
+[:, :, 1] =
+ 1  5
+ 3  7
+
+[:, :, 2] =
+  9  13
+ 11  15
+
+julia> extrema(A, dims = (1,2))
+1×1×2 Array{Tuple{Int64, Int64}, 3}:
+[:, :, 1] =
+ (1, 7)
+
+[:, :, 2] =
+ (9, 15)
+```
+"""
+extrema(A::AbstractArray; dims)
+
+"""
+    extrema(f, A::AbstractArray; dims) -> Array{Tuple}
+
+Compute the minimum and maximum of `f` applied to each element in the given dimensions
+of `A`.
+
+!!! compat "Julia 1.2"
+    This method requires Julia 1.2 or later.
+"""
+extrema(f, A::AbstractArray; dims)
+
+"""
+    extrema!(r, A)
+
+Compute the minimum and maximum value of `A` over the singleton dimensions of `r`, and write results to `r`.
+
+!!! compat "Julia 1.8"
+    This method requires Julia 1.8 or later.
+
+# Examples
+```jldoctest
+julia> A = [1 2; 3 4]
+2×2 Matrix{Int64}:
+ 1  2
+ 3  4
+
+julia> extrema!([(1, 1); (1, 1)], A)
+2-element Vector{Tuple{Int64, Int64}}:
+ (1, 2)
+ (3, 4)
+
+julia> extrema!([(1, 1);; (1, 1)], A)
+1×2 Matrix{Tuple{Int64, Int64}}:
+ (1, 3)  (2, 4)
+```
+"""
+extrema!(r, A)
+
 """
     all(A; dims)
 
@@ -883,7 +986,9 @@ julia> any!([1 1], A)
 any!(r, A)
 
 for (fname, _fname, op) in [(:sum,     :_sum,     :add_sum), (:prod,    :_prod,    :mul_prod),
-                            (:maximum, :_maximum, :max),     (:minimum, :_minimum, :min)]
+                            (:maximum, :_maximum, :max),     (:minimum, :_minimum, :min),
+                            (:extrema, :_extrema, :_extrema_rf)]
+    mapf = fname === :extrema ? :(ExtremaMap(f)) : :f
     @eval begin
         # User-facing methods with keyword arguments
         @inline ($fname)(a::AbstractArray; dims=:, kw...) = ($_fname)(a, dims; kw...)
@@ -891,7 +996,7 @@ for (fname, _fname, op) in [(:sum,     :_sum,     :add_sum), (:prod,    :_prod,
 
         # Underlying implementations using dispatch
         ($_fname)(a, ::Colon; kw...) = ($_fname)(identity, a, :; kw...)
-        ($_fname)(f, a, ::Colon; kw...) = mapreduce(f, $op, a; kw...)
+        ($_fname)(f, a, ::Colon; kw...) = mapreduce($mapf, $op, a; kw...)
     end
 end
 
@@ -904,16 +1009,18 @@ _all(a, ::Colon)                           = _all(identity, a, :)
 
 for (fname, op) in [(:sum, :add_sum), (:prod, :mul_prod),
                     (:maximum, :max), (:minimum, :min),
-                    (:all, :&),       (:any, :|)]
+                    (:all, :&),       (:any, :|),
+                    (:extrema, :_extrema_rf)]
     fname! = Symbol(fname, '!')
     _fname = Symbol('_', fname)
+    mapf = fname === :extrema ? :(ExtremaMap(f)) : :f
     @eval begin
         $(fname!)(f::Function, r::AbstractArray, A::AbstractArray; init::Bool=true) =
-            mapreducedim!(f, $(op), initarray!(r, $(op), init, A), A)
+            mapreducedim!($mapf, $(op), initarray!(r, $mapf, $(op), init, A), A)
         $(fname!)(r::AbstractArray, A::AbstractArray; init::Bool=true) = $(fname!)(identity, r, A; init=init)
 
         $(_fname)(A, dims; kw...)    = $(_fname)(identity, A, dims; kw...)
-        $(_fname)(f, A, dims; kw...) = mapreduce(f, $(op), A; dims=dims, kw...)
+        $(_fname)(f, A, dims; kw...) = mapreduce($mapf, $(op), A; dims=dims, kw...)
     end
 end
 
diff --git a/doc/src/base/collections.md b/doc/src/base/collections.md
index d329ce6ef6119..221adad3954e6 100644
--- a/doc/src/base/collections.md
+++ b/doc/src/base/collections.md
@@ -102,6 +102,7 @@ Base.maximum!
 Base.minimum
 Base.minimum!
 Base.extrema
+Base.extrema!
 Base.argmax
 Base.argmin
 Base.findmax
diff --git a/test/reduce.jl b/test/reduce.jl
index 78ac22c13f366..0e1568b0af901 100644
--- a/test/reduce.jl
+++ b/test/reduce.jl
@@ -244,9 +244,15 @@ prod2(itr) = invoke(prod, Tuple{Any}, itr)
 
 @test_throws "reducing over an empty" maximum(Int[])
 @test_throws "reducing over an empty" minimum(Int[])
+@test_throws "reducing over an empty" extrema(Int[])
 
 @test maximum(Int[]; init=-1) == -1
 @test minimum(Int[]; init=-1) == -1
+@test extrema(Int[]; init=(1, -1)) == (1, -1)
+
+@test maximum(sin, []; init=-1) == -1
+@test minimum(sin, []; init=1) == 1
+@test extrema(sin, []; init=(1, -1)) == (1, -1)
 
 @test maximum(5) == 5
 @test minimum(5) == 5
@@ -257,6 +263,7 @@ let x = [4,3,5,2]
     @test maximum(x) == 5
     @test minimum(x) == 2
     @test extrema(x) == (2, 5)
+    @test Core.Compiler.extrema(x) == (2, 5)
 
     @test maximum(abs2, x) == 25
     @test minimum(abs2, x) == 4
@@ -388,6 +395,9 @@ A = circshift(reshape(1:24,2,3,4), (0,1,1))
         @test maximum(x) === minimum(x) === missing
         @test extrema(x) === (missing, missing)
     end
+    # inputs containing both missing and NaN
+    minimum([NaN;zeros(255);missing]) === missing
+    maximum([NaN;zeros(255);missing]) === missing
 end
 
 # findmin, findmax, argmin, argmax
diff --git a/test/reducedim.jl b/test/reducedim.jl
index db60d353ae95e..512c94d1e2f02 100644
--- a/test/reducedim.jl
+++ b/test/reducedim.jl
@@ -93,6 +93,10 @@ A = Array{Int}(undef, 0, 3)
 @test_throws "reducing over an empty collection is not allowed" maximum(A; dims=1)
 @test maximum(A; dims=1, init=-1) == reshape([-1,-1,-1], 1, 3)
 
+@test maximum(zeros(0, 2); dims=1, init=-1) == fill(-1, 1, 2)
+@test minimum(zeros(0, 2); dims=1, init=1) == ones(1, 2)
+@test extrema(zeros(0, 2); dims=1, init=(1, -1)) == fill((1, -1), 1, 2)
+
 # Test reduction along first dimension; this is special-cased for
 # size(A, 1) >= 16
 Breduc = rand(64, 3)
@@ -454,18 +458,51 @@ end
     for R in (fill(0, 4), fill(0, 4, 1), fill(0, 4, 1, 1))
         @test @inferred(maximum!(R, B)) == reshape(21:24, size(R))
         @test @inferred(minimum!(R, B)) == reshape(1:4, size(R))
+        @test @inferred(extrema!(fill((0,0), size(R)), B)) == reshape(tuple.(1:4, 21:24), size(R))
     end
     for R in (fill(0, 1, 3), fill(0, 1, 3, 1))
         @test @inferred(maximum!(R, B)) == reshape(16:4:24, size(R))
         @test @inferred(minimum!(R, B)) == reshape(1:4:9, size(R))
+        @test @inferred(extrema!(fill((0,0), size(R)), B)) == reshape(tuple.(1:4:9, 16:4:24), size(R))
+    end
+    for (ini, f!) in zip((0,0,(0,0)), (maximum!, minimum!, extrema!))
+        @test_throws DimensionMismatch f!(fill(ini, 4, 1, 1, 1), B)
+        @test_throws DimensionMismatch f!(fill(ini, 1, 3, 1, 1), B)
+        @test_throws DimensionMismatch f!(fill(ini, 1, 1, 2, 1), B)
+    end
+end
+
+function unordered_test_for_extrema(a; dims_test = ((), 1, 2, (1,2), 3))
+    for dims in dims_test
+        vext = extrema(a; dims)
+        vmin, vmax = minimum(a; dims), maximum(a; dims)
+        @test isequal(extrema!(copy(vext), a), vext)
+        @test all(x -> isequal(x[1], x[2:3]), zip(vext,vmin,vmax))
     end
-    @test_throws DimensionMismatch maximum!(fill(0, 4, 1, 1, 1), B)
-    @test_throws DimensionMismatch minimum!(fill(0, 4, 1, 1, 1), B)
-    @test_throws DimensionMismatch maximum!(fill(0, 1, 3, 1, 1), B)
-    @test_throws DimensionMismatch minimum!(fill(0, 1, 3, 1, 1), B)
-    @test_throws DimensionMismatch maximum!(fill(0, 1, 1, 2, 1), B)
-    @test_throws DimensionMismatch minimum!(fill(0, 1, 1, 2, 1), B)
 end
+@testset "0.0,-0.0 test for extrema with dims" begin
+    @test extrema([-0.0;0.0], dims = 1)[1] === (-0.0,0.0)
+    @test tuple(extrema([-0.0;0.0], dims = 2)...) === ((-0.0, -0.0), (0.0, 0.0))
+end
+@testset "NaN/missing test for extrema with dims #43599" begin
+    for sz = (3, 10, 100)
+        for T in (Int, Float64, BigFloat)
+            Aₘ = Matrix{Union{T, Missing}}(rand(-sz:sz, sz, sz))
+            Aₘ[rand(1:sz*sz, sz)] .= missing
+            unordered_test_for_extrema(Aₘ)
+            if T <: AbstractFloat
+                Aₙ = map(i -> ismissing(i) ? T(NaN) : i, Aₘ)
+                unordered_test_for_extrema(Aₙ)
+                p = rand(1:sz*sz, sz)
+                Aₘ[p] .= NaN
+                unordered_test_for_extrema(Aₘ)
+            end
+        end
+    end
+end
+@test_broken minimum([missing;BigInt(1)], dims = 1)
+@test_broken maximum([missing;BigInt(1)], dims = 1)
+@test_broken extrema([missing;BigInt(1)], dims = 1)
 
 # issue #26709
 @testset "dimensional reduce with custom non-bitstype types" begin
@@ -510,3 +547,11 @@ end
 @testset "type stability (issue #43461)" begin
     @test (@inferred maximum(Float64, reshape(1:4,2,:); dims = 2)) == reshape([3,4],2,1)
 end
+
+@testset "Min/Max initialization test" begin
+    A = Vector{Union{Missing,Int}}(1:4)
+    A[2] = missing
+    @test_broken @inferred(minimum(exp, A; dims = 1))[1] === missing
+    @test_broken @inferred(maximum(exp, A; dims = 1))[1] === missing
+    @test_broken @inferred(extrema(exp, A; dims = 1))[1] === (missing, missing)
+end