From 7b0e72efc8f7fbf2c0bd5d30a8b6e2e22288547f Mon Sep 17 00:00:00 2001
From: Takafumi Arakaki <aka.tkf@gmail.com>
Date: Wed, 10 Jun 2020 23:15:39 -0700
Subject: [PATCH 01/19] Define `extrema` using `mapreduce`; support `init`

---
 base/compiler/compiler.jl | 18 ++++++++++
 base/multidimensional.jl  | 37 ++++++++++----------
 base/operators.jl         | 52 ----------------------------
 base/reduce.jl            | 71 ++++++++++++++++++++++++++++++++++++++-
 test/reduce.jl            |  6 ++++
 test/reducedim.jl         |  4 +++
 6 files changed, 117 insertions(+), 71 deletions(-)

diff --git a/base/compiler/compiler.jl b/base/compiler/compiler.jl
index 4ac2e3ca9a4fd..0a2f2ba0e4895 100644
--- a/base/compiler/compiler.jl
+++ b/base/compiler/compiler.jl
@@ -48,6 +48,24 @@ include("operators.jl")
 include("pointer.jl")
 include("refvalue.jl")
 
+# required for bootstrap
+extrema(itr) = extrema(identity, itr)
+function extrema(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
+
 # checked arithmetic
 const checked_add = +
 const checked_sub = -
diff --git a/base/multidimensional.jl b/base/multidimensional.jl
index 51de7184cb92f..ec68f55a3a864 100644
--- a/base/multidimensional.jl
+++ b/base/multidimensional.jl
@@ -1566,9 +1566,13 @@ _unique_dims(A::AbstractArray, dims::Colon) = invoke(unique, Tuple{Any}, A)
 end
 
 """
-    extrema(A::AbstractArray; dims) -> Array{Tuple}
+    extrema([f,] A::AbstractArray; dims, [init]) -> Array{Tuple}
 
-Compute the minimum and maximum elements of an array over the given dimensions.
+Compute the minimum and maximum of `f` applied to each element (if given) in the given
+dimensions of `A`.
+
+!!! compat "Julia 1.2"
+    2-argument `extrema` method requires Julia 1.2 or later.
 
 # Examples
 ```jldoctest
@@ -1591,22 +1595,19 @@ julia> extrema(A, dims = (1,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)
+extrema(f::F, A::AbstractArray; dims=:, init=_InitialValue()) where {F} =
+    _extrema_dims(f, A, dims, init)
+
+_extrema_dims(f::F, A::AbstractArray, ::Colon, init) where {F} =
+    mapreduce(_DupY(f), _extrema_rf, A; init = init)
+_extrema_dims(f::F, A::AbstractArray, ::Colon, ::_InitialValue) where {F} =
+    mapreduce(_DupY(f), _extrema_rf, A)
+# Note: not passing `init = _InitialValue()` since user-defined
+# `reduce`/`foldl` could cannot be aware of `Base._InitialValue`.
+
+_extrema_dims(f::F, A::AbstractArray, dims, init) where {F} =
+    mapreduce(_DupY(f), _extrema_rf, A; dims = dims, init = init)
+function _extrema_dims(f::F, A::AbstractArray, dims, ::_InitialValue) where {F}
     sz = [size(A)...]
     for d in dims
         sz[d] = 1
diff --git a/base/operators.jl b/base/operators.jl
index 605c0775ad00f..c59c3a18d6377 100644
--- a/base/operators.jl
+++ b/base/operators.jl
@@ -443,58 +443,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 c73caaeac7a97..b850ab9c7577a 100644
--- a/base/reduce.jl
+++ b/base/reduce.jl
@@ -581,7 +581,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)
@@ -762,6 +762,75 @@ Inf
 """
 minimum(a; kw...) = mapreduce(identity, min, a; kw...)
 
+"""
+    extrema(itr; [init]) -> Tuple
+
+Compute both the minimum and maximum 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/second element is a neutral element for `min`/`max` (i.e. which is greater/less than
+or equal to any other element). This is required because it is unspecified whether `init`
+is used for non-empty collections. Note: it implies that, for empty `itr`, the first
+element is typically _greater_ than the last element. This is a "paradoxical" but yet
+expected result.
+
+!!! compat "Julia 1.6"
+    Keyword argument `init` requires Julia 1.6 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]) -> 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`.
+
+The value returned for empty `itr` can be specified by `init`. It must be a 2-tuple whose
+first/second element is a neutral element for `min`/`max` (i.e. which is greater/less than
+or equal to any other element). This is required because it is unspecified whether `init`
+is used for non-empty collections. Note: it implies that, for empty `itr`, the first
+element is typically _greater_ than the last element. This is a "paradoxical" but yet
+expected result.
+
+!!! compat "Julia 1.2"
+    This method requires Julia 1.2 or later.
+
+!!! compat "Julia 1.6"
+    Keyword argument `init` requires Julia 1.6 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(_DupY(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 _DupY{F} <: Function
+    f::F
+end
+_DupY(f::Type{T}) where {T} = _DupY{Type{T}}(f)
+@inline (f::_DupY)(x) = (y = f.f(x); (y, y))
+
+@inline _extrema_rf((min1, max1), (min2, max2)) = (min(min1, min2), max(max1, max2))
+
 ## all & any
 
 """
diff --git a/test/reduce.jl b/test/reduce.jl
index 4688709b099c9..6f39fd9f6d3c4 100644
--- a/test/reduce.jl
+++ b/test/reduce.jl
@@ -241,9 +241,15 @@ prod2(itr) = invoke(prod, Tuple{Any}, itr)
 
 @test_throws ArgumentError maximum(Int[])
 @test_throws ArgumentError minimum(Int[])
+@test_throws ArgumentError 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
diff --git a/test/reducedim.jl b/test/reducedim.jl
index e51a075496f1c..b0668745d55a6 100644
--- a/test/reducedim.jl
+++ b/test/reducedim.jl
@@ -80,6 +80,10 @@ A = Array{Int}(undef, 0, 3)
 @test_throws ArgumentError 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)

From c6ffc8739f572bf3fc57eb404474f727509c76e4 Mon Sep 17 00:00:00 2001
From: Takafumi Arakaki <aka.tkf@gmail.com>
Date: Fri, 12 Jun 2020 21:19:00 -0700
Subject: [PATCH 02/19] Add NEWS

---
 NEWS.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/NEWS.md b/NEWS.md
index b6ffe8c872af5..12d0609882366 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -47,7 +47,7 @@ Standard library changes
 * The function `isapprox(x,y)` now accepts the `norm` keyword argument also for numeric (i.e., non-array) arguments `x` and `y` ([#35883]).
 * `view`, `@view`, and `@views` now work on `AbstractString`s, returning a `SubString` when appropriate ([#35879]).
 * All `AbstractUnitRange{<:Integer}`s now work with `SubString`, `view`, `@view` and `@views` on strings ([#35879]).
-* `sum`, `prod`, `maximum`, and `minimum` now support `init` keyword argument ([#36188], [#35839]).
+* `sum`, `prod`, `maximum`, `minimum`, and `extrema` now support `init` keyword argument ([#36188], [#35839], [#36265]).
 
 #### LinearAlgebra
 * New method `LinearAlgebra.issuccess(::CholeskyPivoted)` for checking whether pivoted Cholesky factorization was successful ([#36002]).

From 194c329435e7a312bfb22bcaa3d852c76f111c20 Mon Sep 17 00:00:00 2001
From: Takafumi Arakaki <aka.tkf@gmail.com>
Date: Fri, 12 Jun 2020 21:48:53 -0700
Subject: [PATCH 03/19] Fix a typo

---
 base/multidimensional.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/base/multidimensional.jl b/base/multidimensional.jl
index ec68f55a3a864..97482fb5124e1 100644
--- a/base/multidimensional.jl
+++ b/base/multidimensional.jl
@@ -1603,7 +1603,7 @@ _extrema_dims(f::F, A::AbstractArray, ::Colon, init) where {F} =
 _extrema_dims(f::F, A::AbstractArray, ::Colon, ::_InitialValue) where {F} =
     mapreduce(_DupY(f), _extrema_rf, A)
 # Note: not passing `init = _InitialValue()` since user-defined
-# `reduce`/`foldl` could cannot be aware of `Base._InitialValue`.
+# `reduce`/`foldl` cannot be aware of `Base._InitialValue`.
 
 _extrema_dims(f::F, A::AbstractArray, dims, init) where {F} =
     mapreduce(_DupY(f), _extrema_rf, A; dims = dims, init = init)

From 793a02056357edb590e01589a6aaab8ca5cd82b4 Mon Sep 17 00:00:00 2001
From: Takafumi Arakaki <aka.tkf@gmail.com>
Date: Fri, 12 Jun 2020 22:42:15 -0700
Subject: [PATCH 04/19] Reword a comment slightly

---
 base/multidimensional.jl | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/base/multidimensional.jl b/base/multidimensional.jl
index 97482fb5124e1..119f6d7ca805f 100644
--- a/base/multidimensional.jl
+++ b/base/multidimensional.jl
@@ -1603,7 +1603,8 @@ _extrema_dims(f::F, A::AbstractArray, ::Colon, init) where {F} =
 _extrema_dims(f::F, A::AbstractArray, ::Colon, ::_InitialValue) where {F} =
     mapreduce(_DupY(f), _extrema_rf, A)
 # Note: not passing `init = _InitialValue()` since user-defined
-# `reduce`/`foldl` cannot be aware of `Base._InitialValue`.
+# `reduce`/`foldl` cannot be aware of `Base._InitialValue` that is an
+# internal implementation detail.
 
 _extrema_dims(f::F, A::AbstractArray, dims, init) where {F} =
     mapreduce(_DupY(f), _extrema_rf, A; dims = dims, init = init)

From 5ebf8ad4ec97c2dea49401373751782438eeed45 Mon Sep 17 00:00:00 2001
From: Takafumi Arakaki <aka.tkf@gmail.com>
Date: Fri, 19 Jun 2020 12:44:15 -0700
Subject: [PATCH 05/19] Apply suggestions from code review

Co-authored-by: Milan Bouchet-Valat <nalimilan@club.fr>
---
 base/multidimensional.jl | 6 +++---
 base/reduce.jl           | 8 ++++----
 2 files changed, 7 insertions(+), 7 deletions(-)

diff --git a/base/multidimensional.jl b/base/multidimensional.jl
index 119f6d7ca805f..aa4c0526e6581 100644
--- a/base/multidimensional.jl
+++ b/base/multidimensional.jl
@@ -1568,11 +1568,11 @@ end
 """
     extrema([f,] A::AbstractArray; dims, [init]) -> Array{Tuple}
 
-Compute the minimum and maximum of `f` applied to each element (if given) in the given
-dimensions of `A`.
+Compute the minimum and maximum elements of `A` over dimensions `dims`.
+If `f` is provided, return the minimum and maximum elements after applying `f` to them.
 
 !!! compat "Julia 1.2"
-    2-argument `extrema` method requires Julia 1.2 or later.
+    The `extrema(f, A)` method requires Julia 1.2 or later.
 
 # Examples
 ```jldoctest
diff --git a/base/reduce.jl b/base/reduce.jl
index b850ab9c7577a..44032ce8523af 100644
--- a/base/reduce.jl
+++ b/base/reduce.jl
@@ -768,8 +768,8 @@ minimum(a; kw...) = mapreduce(identity, min, a; kw...)
 Compute both the minimum and maximum 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/second element is a neutral element for `min`/`max` (i.e. which is greater/less than
-or equal to any other element). This is required because it is unspecified whether `init`
+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).
 is used for non-empty collections. Note: it implies that, for empty `itr`, the first
 element is typically _greater_ than the last element. This is a "paradoxical" but yet
 expected result.
@@ -798,8 +798,8 @@ Compute both the minimum and maximum of `f` applied to each element in `itr` and
 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/second element is a neutral element for `min`/`max` (i.e. which is greater/less than
-or equal to any other element). This is required because it is unspecified whether `init`
+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).
 is used for non-empty collections. Note: it implies that, for empty `itr`, the first
 element is typically _greater_ than the last element. This is a "paradoxical" but yet
 expected result.

From 2d70b7d5300e931a301ad0c828d01f7e781e768e Mon Sep 17 00:00:00 2001
From: Takafumi Arakaki <aka.tkf@gmail.com>
Date: Fri, 19 Jun 2020 12:46:48 -0700
Subject: [PATCH 06/19] Fixup docstrings

---
 base/reduce.jl | 14 ++++++--------
 1 file changed, 6 insertions(+), 8 deletions(-)

diff --git a/base/reduce.jl b/base/reduce.jl
index 44032ce8523af..15b0715d55aa2 100644
--- a/base/reduce.jl
+++ b/base/reduce.jl
@@ -769,10 +769,9 @@ Compute both the minimum and maximum element in a single pass, and return them a
 
 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).
-is used for non-empty collections. Note: it implies that, for empty `itr`, the first
-element is typically _greater_ than the last element. This is a "paradoxical" but yet
-expected result.
+(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 first element is typically
+_greater_ than the last element. This is a "paradoxical" but yet expected result.
 
 !!! compat "Julia 1.6"
     Keyword argument `init` requires Julia 1.6 or later.
@@ -799,10 +798,9 @@ 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).
-is used for non-empty collections. Note: it implies that, for empty `itr`, the first
-element is typically _greater_ than the last element. This is a "paradoxical" but yet
-expected result.
+(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 first element is typically
+_greater_ than the last element. This is a "paradoxical" but yet expected result.
 
 !!! compat "Julia 1.2"
     This method requires Julia 1.2 or later.

From 317548a457f5d91e9602dfdcb8f0b039e515d8d5 Mon Sep 17 00:00:00 2001
From: Takafumi Arakaki <aka.tkf@gmail.com>
Date: Fri, 19 Jun 2020 13:30:44 -0700
Subject: [PATCH 07/19] Add simple tests for Core.Compiler.extrema

---
 test/reduce.jl | 3 +++
 1 file changed, 3 insertions(+)

diff --git a/test/reduce.jl b/test/reduce.jl
index 6f39fd9f6d3c4..2ea40f5289c9a 100644
--- a/test/reduce.jl
+++ b/test/reduce.jl
@@ -255,15 +255,18 @@ prod2(itr) = invoke(prod, Tuple{Any}, itr)
 @test minimum(5) == 5
 @test extrema(5) == (5, 5)
 @test extrema(abs2, 5) == (25, 25)
+@test Core.Compiler.extrema(abs2, 5) == (25, 25)
 
 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
     @test extrema(abs2, x) == (4, 25)
+    @test Core.Compiler.extrema(abs2, x) == (4, 25)
 end
 
 @test maximum([-0.,0.]) === 0.0

From 360b3463cef1f03681ffdbffbb6cb005f9ac09f4 Mon Sep 17 00:00:00 2001
From: Takafumi Arakaki <aka.tkf@gmail.com>
Date: Sat, 20 Jun 2020 01:19:54 -0700
Subject: [PATCH 08/19] Name the output tuple elements (mn, mx)

---
 base/reduce.jl | 21 ++++++++++++---------
 1 file changed, 12 insertions(+), 9 deletions(-)

diff --git a/base/reduce.jl b/base/reduce.jl
index 15b0715d55aa2..335de05f94d6d 100644
--- a/base/reduce.jl
+++ b/base/reduce.jl
@@ -763,15 +763,17 @@ Inf
 minimum(a; kw...) = mapreduce(identity, min, a; kw...)
 
 """
-    extrema(itr; [init]) -> Tuple
+    extrema(itr; [init]) -> (mn, mx)
 
-Compute both the minimum and maximum element in a single pass, and return them as a 2-tuple.
+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). It is used for non-empty
-collections. Note: it implies that, for empty `itr`, the first element is typically
-_greater_ than the last element. This is a "paradoxical" but yet expected result.
+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.6"
     Keyword argument `init` requires Julia 1.6 or later.
@@ -791,16 +793,17 @@ julia> extrema([]; init = (Inf, -Inf))
 extrema(itr; kw...) = extrema(identity, itr; kw...)
 
 """
-    extrema(f, itr; [init]) -> Tuple
+    extrema(f, itr; [init]) -> (mn, mx)
 
-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`.
+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 first element is typically
-_greater_ than the last element. This is a "paradoxical" but yet expected result.
+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.

From f27f874ed0babf750ff9dd0955f52c49a2ea8618 Mon Sep 17 00:00:00 2001
From: Takafumi Arakaki <aka.tkf@gmail.com>
Date: Fri, 12 Nov 2021 15:44:19 -0800
Subject: [PATCH 09/19] Update base/multidimensional.jl

---
 base/multidimensional.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/base/multidimensional.jl b/base/multidimensional.jl
index abdfd16c23bc4..c22d6cd922529 100644
--- a/base/multidimensional.jl
+++ b/base/multidimensional.jl
@@ -1752,7 +1752,7 @@ _extrema_dims(f::F, A::AbstractArray, ::Colon, ::_InitialValue) where {F} =
 _extrema_dims(f::F, A::AbstractArray, dims, init) where {F} =
     mapreduce(_DupY(f), _extrema_rf, A; dims = dims, init = init)
 function _extrema_dims(f::F, A::AbstractArray, dims, ::_InitialValue) where {F}
-    sz = [size(A)...]
+    sz = size(A)
     for d in dims
         sz = setindex(sz, 1, d)
     end

From ce07cbf185295bd21dedac70ccabeeb3b4fcf22d Mon Sep 17 00:00:00 2001
From: Takafumi Arakaki <aka.tkf@gmail.com>
Date: Fri, 12 Nov 2021 18:54:56 -0500
Subject: [PATCH 10/19] Fix compat annotations

---
 base/reduce.jl | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/base/reduce.jl b/base/reduce.jl
index 1106d5409348e..832ae35bf04c9 100644
--- a/base/reduce.jl
+++ b/base/reduce.jl
@@ -798,8 +798,8 @@ collections. Note: it implies that, for empty `itr`, the returned value `(mn, mx
 `mn ≥ mx` even though for non-empty `itr` it  satisfies `mn ≤ mx`.  This is a "paradoxical"
 but yet expected result.
 
-!!! compat "Julia 1.6"
-    Keyword argument `init` requires Julia 1.6 or later.
+!!! compat "Julia 1.8"
+    Keyword argument `init` requires Julia 1.8 or later.
 
 # Examples
 ```jldoctest
@@ -831,8 +831,8 @@ but yet expected result.
 !!! compat "Julia 1.2"
     This method requires Julia 1.2 or later.
 
-!!! compat "Julia 1.6"
-    Keyword argument `init` requires Julia 1.6 or later.
+!!! compat "Julia 1.8"
+    Keyword argument `init` requires Julia 1.8 or later.
 
 # Examples
 ```jldoctest

From 3ceaaae94bb5a26c36a72c5e824f90c1821ca8a4 Mon Sep 17 00:00:00 2001
From: Takafumi Arakaki <aka.tkf@gmail.com>
Date: Fri, 12 Nov 2021 18:59:07 -0500
Subject: [PATCH 11/19] Improve extrema docstring

Co-authored-by: Tim Holy <tim.holy@gmail.com>
---
 base/reduce.jl | 7 +++----
 1 file changed, 3 insertions(+), 4 deletions(-)

diff --git a/base/reduce.jl b/base/reduce.jl
index 832ae35bf04c9..23a56c0c5f45f 100644
--- a/base/reduce.jl
+++ b/base/reduce.jl
@@ -793,10 +793,9 @@ 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). 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.
+(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.

From 5669d6349569ff83394208609c1d5c4b2f86fa17 Mon Sep 17 00:00:00 2001
From: Takafumi Arakaki <aka.tkf@gmail.com>
Date: Fri, 12 Nov 2021 23:34:05 -0500
Subject: [PATCH 12/19] Fix tests for SparseArrays

---
 stdlib/SparseArrays/test/higherorderfns.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/stdlib/SparseArrays/test/higherorderfns.jl b/stdlib/SparseArrays/test/higherorderfns.jl
index 8da605cf6c0c0..59079682938ee 100644
--- a/stdlib/SparseArrays/test/higherorderfns.jl
+++ b/stdlib/SparseArrays/test/higherorderfns.jl
@@ -709,8 +709,8 @@ end
     @test extrema(f, x) == extrema(f, y)
     @test extrema(spzeros(n, n)) == (0.0, 0.0)
     @test extrema(spzeros(n)) == (0.0, 0.0)
-    @test_throws ArgumentError extrema(spzeros(0, 0))
-    @test_throws ArgumentError extrema(spzeros(0))
+    @test_throws "reducing over an empty" extrema(spzeros(0, 0))
+    @test_throws "reducing over an empty" extrema(spzeros(0))
     @test extrema(sparse(ones(n, n))) == (1.0, 1.0)
     @test extrema(sparse(ones(n))) == (1.0, 1.0)
     @test extrema(A; dims=:) == extrema(B; dims=:)

From ce163962219299d47d24dcaddba4489ff30d55da Mon Sep 17 00:00:00 2001
From: N5N3 <2642243996@qq.com>
Date: Sat, 15 Jan 2022 12:14:11 +0800
Subject: [PATCH 13/19] API clean

and export `extrema!`

fix doc

Update collections.md

drop sparsearray change

Co-Authored-By: Takafumi Arakaki <29282+tkf@users.noreply.github.com>
---
 base/compiler/compiler.jl                  | 31 +++----
 base/exports.jl                            |  1 +
 base/mpfr.jl                               | 11 +++
 base/multidimensional.jl                   | 76 -----------------
 base/reduce.jl                             | 16 +++-
 base/reducedim.jl                          | 94 +++++++++++++++++++---
 doc/src/base/collections.md                |  1 +
 stdlib/SparseArrays/test/higherorderfns.jl |  4 +-
 test/reduce.jl                             |  2 -
 9 files changed, 123 insertions(+), 113 deletions(-)

diff --git a/base/compiler/compiler.jl b/base/compiler/compiler.jl
index 6a8ea79c2142a..aafd80c1127ff 100644
--- a/base/compiler/compiler.jl
+++ b/base/compiler/compiler.jl
@@ -55,24 +55,6 @@ include("operators.jl")
 include("pointer.jl")
 include("refvalue.jl")
 
-# required for bootstrap
-extrema(itr) = extrema(identity, itr)
-function extrema(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
-
 # checked arithmetic
 const checked_add = +
 const checked_sub = -
@@ -149,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 8174488897cb5..03184bceed68b 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/mpfr.jl b/base/mpfr.jl
index e85f281619ac0..98158e6f44121 100644
--- a/base/mpfr.jl
+++ b/base/mpfr.jl
@@ -706,6 +706,17 @@ function min(x::BigFloat, y::BigFloat)
     return z
 end
 
+# We don't want allocations during `extrema`.
+# TODO: avoid allocations in `min/max`
+function Base._extrema_rf(x::NTuple{2,BigFloat}, y::NTuple{2,BigFloat})
+    (x1, x2), (y1, y2) = x, y
+    isnan(x1) && return x
+    isnan(y1) && return y
+    z1 = x1 < y1 || signbit(x1) > signbit(y1) ? x1 : y1
+    z2 = x2 < y2 || signbit(x2) > signbit(y2) ? y2 : x2
+    z1, z2
+end
+
 function modf(x::BigFloat)
     zint = BigFloat()
     zfloat = BigFloat()
diff --git a/base/multidimensional.jl b/base/multidimensional.jl
index 4a1a86512c769..7a504e19b63d0 100644
--- a/base/multidimensional.jl
+++ b/base/multidimensional.jl
@@ -1702,82 +1702,6 @@ _unique_dims(A::AbstractArray, dims::Colon) = invoke(unique, Tuple{Any}, A)
     end
 end
 
-"""
-    extrema([f,] A::AbstractArray; dims, [init]) -> Array{Tuple}
-
-Compute the minimum and maximum elements of `A` over dimensions `dims`.
-If `f` is provided, return the minimum and maximum elements after applying `f` to them.
-
-!!! compat "Julia 1.2"
-    The `extrema(f, A)` method requires Julia 1.2 or later.
-
-# 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(f::F, A::AbstractArray; dims=:, init=_InitialValue()) where {F} =
-    _extrema_dims(f, A, dims, init)
-
-_extrema_dims(f::F, A::AbstractArray, ::Colon, init) where {F} =
-    mapreduce(_DupY(f), _extrema_rf, A; init = init)
-_extrema_dims(f::F, A::AbstractArray, ::Colon, ::_InitialValue) where {F} =
-    mapreduce(_DupY(f), _extrema_rf, A)
-# Note: not passing `init = _InitialValue()` since user-defined
-# `reduce`/`foldl` cannot be aware of `Base._InitialValue` that is an
-# internal implementation detail.
-
-_extrema_dims(f::F, A::AbstractArray, dims, init) where {F} =
-    mapreduce(_DupY(f), _extrema_rf, A; dims = dims, init = init)
-function _extrema_dims(f::F, A::AbstractArray, dims, ::_InitialValue) where {F}
-    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/reduce.jl b/base/reduce.jl
index 5f02a8917d86f..808e7c348a418 100644
--- a/base/reduce.jl
+++ b/base/reduce.jl
@@ -842,17 +842,25 @@ julia> extrema(sin, Real[]; init = (1.0, -1.0))  # good, since -1 ≤ sin(::Real
 (1.0, -1.0)
 ```
 """
-extrema(f, itr; kw...) = mapreduce(_DupY(f), _extrema_rf, itr; kw...)
+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 _DupY{F} <: Function
+struct ExtremaMap{F} <: Function
     f::F
 end
-_DupY(f::Type{T}) where {T} = _DupY{Type{T}}(f)
-@inline (f::_DupY)(x) = (y = f.f(x); (y, y))
+ExtremaMap(::Type{T}) where {T} = ExtremaMap{Type{T}}(T)
+@inline (f::ExtremaMap)(x) = (y = f.f(x); (y, y))
 
 @inline _extrema_rf((min1, max1), (min2, max2)) = (min(min1, min2), max(max1, max2))
+# optimization for AbstractFloat
+function _extrema_rf(x::NTuple{2,T}, y::NTuple{2,T}) where {T<:AbstractFloat}
+    (x1, x2), (y1, y2) = x, y
+    anynan = isnan(x1)|isnan(y1)
+    z1 = ifelse(anynan, x1-y1, ifelse(signbit(x1-y1), x1, y1))
+    z2 = ifelse(anynan, x1-y1, ifelse(signbit(x2-y2), y2, x2))
+    z1, z2
+end
 
 ## findmax, findmin, argmax & argmin
 
diff --git a/base/reducedim.jl b/base/reducedim.jl
index 02b2345038bcc..565b61ba9e5c3 100644
--- a/base/reducedim.jl
+++ b/base/reducedim.jl
@@ -77,15 +77,15 @@ 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
+# for min/max copyfirst is not correct for initialization, use `mapfirst!` instead
+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
@@ -435,7 +435,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 +737,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 +951,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 +961,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 +974,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/stdlib/SparseArrays/test/higherorderfns.jl b/stdlib/SparseArrays/test/higherorderfns.jl
index 59079682938ee..8da605cf6c0c0 100644
--- a/stdlib/SparseArrays/test/higherorderfns.jl
+++ b/stdlib/SparseArrays/test/higherorderfns.jl
@@ -709,8 +709,8 @@ end
     @test extrema(f, x) == extrema(f, y)
     @test extrema(spzeros(n, n)) == (0.0, 0.0)
     @test extrema(spzeros(n)) == (0.0, 0.0)
-    @test_throws "reducing over an empty" extrema(spzeros(0, 0))
-    @test_throws "reducing over an empty" extrema(spzeros(0))
+    @test_throws ArgumentError extrema(spzeros(0, 0))
+    @test_throws ArgumentError extrema(spzeros(0))
     @test extrema(sparse(ones(n, n))) == (1.0, 1.0)
     @test extrema(sparse(ones(n))) == (1.0, 1.0)
     @test extrema(A; dims=:) == extrema(B; dims=:)
diff --git a/test/reduce.jl b/test/reduce.jl
index 155ce507a0c0d..18adfc9d4e10d 100644
--- a/test/reduce.jl
+++ b/test/reduce.jl
@@ -258,7 +258,6 @@ prod2(itr) = invoke(prod, Tuple{Any}, itr)
 @test minimum(5) == 5
 @test extrema(5) == (5, 5)
 @test extrema(abs2, 5) == (25, 25)
-@test Core.Compiler.extrema(abs2, 5) == (25, 25)
 
 let x = [4,3,5,2]
     @test maximum(x) == 5
@@ -269,7 +268,6 @@ let x = [4,3,5,2]
     @test maximum(abs2, x) == 25
     @test minimum(abs2, x) == 4
     @test extrema(abs2, x) == (4, 25)
-    @test Core.Compiler.extrema(abs2, x) == (4, 25)
 end
 
 @test maximum([-0.,0.]) === 0.0

From cfb2e8b4b0d43e092e60f4911cc529071c2ade61 Mon Sep 17 00:00:00 2001
From: N5N3 <2642243996@qq.com>
Date: Sat, 8 Jan 2022 14:51:16 +0800
Subject: [PATCH 14/19] add test for extrema(!)

---
 test/reducedim.jl | 50 +++++++++++++++++++++++++++++++++++++++++------
 1 file changed, 44 insertions(+), 6 deletions(-)

diff --git a/test/reducedim.jl b/test/reducedim.jl
index c5995b1378237..111733c2c8c27 100644
--- a/test/reducedim.jl
+++ b/test/reducedim.jl
@@ -458,17 +458,47 @@ 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
+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, BigInt, 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
-    @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
 
 # issue #26709
@@ -514,3 +544,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 @inferred(minimum(exp, A; dims = 1))[1] === missing
+    @test @inferred(maximum(exp, A; dims = 1))[1] === missing
+    @test @inferred(extrema(exp, A; dims = 1))[1] === (missing, missing)
+end

From d9fce2a8b3534636fe50deb12ebd533f9df1bcd9 Mon Sep 17 00:00:00 2001
From: N5N3 <2642243996@qq.com>
Date: Sun, 9 Jan 2022 12:03:01 +0800
Subject: [PATCH 15/19] remove short circuit for correctness

This optimization gives wrong result, if `NaN` and `missing` both exist.

add missing test
---
 base/reduce.jl | 5 -----
 test/reduce.jl | 3 +++
 2 files changed, 3 insertions(+), 5 deletions(-)

diff --git a/base/reduce.jl b/base/reduce.jl
index 808e7c348a418..4349ec3088589 100644
--- a/base/reduce.jl
+++ b/base/reduce.jl
@@ -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]))
diff --git a/test/reduce.jl b/test/reduce.jl
index 18adfc9d4e10d..0e1568b0af901 100644
--- a/test/reduce.jl
+++ b/test/reduce.jl
@@ -395,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

From 065a28f739e3cb3a97fc5a1378d81561d17c1065 Mon Sep 17 00:00:00 2001
From: N5N3 <2642243996@qq.com>
Date: Sun, 9 Jan 2022 10:25:45 +0800
Subject: [PATCH 16/19] Remove `_extrema_rf` optimization for `AbstractFloat`
 (follow advise)

Add TODO

apply suggestions

Update NEWS.md

Co-Authored-By: Takafumi Arakaki <takafumi.a@gmail.com>
Co-Authored-By: Takafumi Arakaki <29282+tkf@users.noreply.github.com>
---
 NEWS.md           |  2 +-
 base/mpfr.jl      | 11 -----------
 base/reduce.jl    |  9 +--------
 base/reducedim.jl |  1 -
 4 files changed, 2 insertions(+), 21 deletions(-)

diff --git a/NEWS.md b/NEWS.md
index bdc8b4c61fc84..4e8b34e66589c 100644
--- a/NEWS.md
+++ b/NEWS.md
@@ -96,7 +96,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]).
+* `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/mpfr.jl b/base/mpfr.jl
index 98158e6f44121..e85f281619ac0 100644
--- a/base/mpfr.jl
+++ b/base/mpfr.jl
@@ -706,17 +706,6 @@ function min(x::BigFloat, y::BigFloat)
     return z
 end
 
-# We don't want allocations during `extrema`.
-# TODO: avoid allocations in `min/max`
-function Base._extrema_rf(x::NTuple{2,BigFloat}, y::NTuple{2,BigFloat})
-    (x1, x2), (y1, y2) = x, y
-    isnan(x1) && return x
-    isnan(y1) && return y
-    z1 = x1 < y1 || signbit(x1) > signbit(y1) ? x1 : y1
-    z2 = x2 < y2 || signbit(x2) > signbit(y2) ? y2 : x2
-    z1, z2
-end
-
 function modf(x::BigFloat)
     zint = BigFloat()
     zfloat = BigFloat()
diff --git a/base/reduce.jl b/base/reduce.jl
index 4349ec3088589..13e1b69c79ede 100644
--- a/base/reduce.jl
+++ b/base/reduce.jl
@@ -847,15 +847,8 @@ 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))
-# optimization for AbstractFloat
-function _extrema_rf(x::NTuple{2,T}, y::NTuple{2,T}) where {T<:AbstractFloat}
-    (x1, x2), (y1, y2) = x, y
-    anynan = isnan(x1)|isnan(y1)
-    z1 = ifelse(anynan, x1-y1, ifelse(signbit(x1-y1), x1, y1))
-    z2 = ifelse(anynan, x1-y1, ifelse(signbit(x2-y2), y2, x2))
-    z1, z2
-end
 
 ## findmax, findmin, argmax & argmin
 
diff --git a/base/reducedim.jl b/base/reducedim.jl
index 565b61ba9e5c3..1118cb56a678e 100644
--- a/base/reducedim.jl
+++ b/base/reducedim.jl
@@ -80,7 +80,6 @@ for (Op, initfun) in ((:(typeof(add_sum)), :zero), (:(typeof(mul_prod)), :one))
     @eval initarray!(a::AbstractArray{T}, ::Any, ::$(Op), init::Bool, src::AbstractArray) where {T} = (init && fill!(a, $(initfun)(T)); a)
 end
 
-# for min/max copyfirst is not correct for initialization, use `mapfirst!` instead
 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)
 

From b6a93bdd1b84a5237fe295fb9d50ec19e24ee0e9 Mon Sep 17 00:00:00 2001
From: N5N3 <2642243996@qq.com>
Date: Sun, 16 Jan 2022 14:53:08 +0800
Subject: [PATCH 17/19] Re-implement `reducedim_init` for extrema

---
 base/reducedim.jl | 35 +++++++++++++++++++++++++++++++++++
 1 file changed, 35 insertions(+)

diff --git a/base/reducedim.jl b/base/reducedim.jl
index 1118cb56a678e..1b437c556d4f3 100644
--- a/base/reducedim.jl
+++ b/base/reducedim.jl
@@ -164,6 +164,41 @@ 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)))
+    Tr = v0[1] isa T && v0[2] isa T ? NTuple{2,T} : typeof(v0)
+
+    # 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
+        T1, T2 = Tr.parameters
+        # TODO: Some types, like BigInt, don't support typemin/typemax.
+        # So a Matrix{Union{BigInt, Missing}} can still error here.
+        v0 = typemax(nonmissingtype(T1)), typemin(nonmissingtype(T2))
+    end
+    # v0 may have changed type.
+    Tr = v0[1] isa T && v0[2] isa T ? NTuple{2,T} : typeof(v0)
+
+    return reducedim_initarray(A, region, v0, Tr)
+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))
 

From a2a039d2ad907783f842d72ce5acd39c32c3ad5d Mon Sep 17 00:00:00 2001
From: N5N3 <2642243996@qq.com>
Date: Sun, 16 Jan 2022 15:32:16 +0800
Subject: [PATCH 18/19] Test fix

Mark `BigInt` as broken
---
 test/reducedim.jl | 11 +++++++----
 1 file changed, 7 insertions(+), 4 deletions(-)

diff --git a/test/reducedim.jl b/test/reducedim.jl
index 111733c2c8c27..512c94d1e2f02 100644
--- a/test/reducedim.jl
+++ b/test/reducedim.jl
@@ -486,7 +486,7 @@ end
 end
 @testset "NaN/missing test for extrema with dims #43599" begin
     for sz = (3, 10, 100)
-        for T in (Int, BigInt, Float64, BigFloat)
+        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ₘ)
@@ -500,6 +500,9 @@ 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
@@ -548,7 +551,7 @@ end
 @testset "Min/Max initialization test" begin
     A = Vector{Union{Missing,Int}}(1:4)
     A[2] = missing
-    @test @inferred(minimum(exp, A; dims = 1))[1] === missing
-    @test @inferred(maximum(exp, A; dims = 1))[1] === missing
-    @test @inferred(extrema(exp, A; dims = 1))[1] === (missing, 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

From 98d27434ca619c3b1bd9c2dd2cf7786550c38efc Mon Sep 17 00:00:00 2001
From: N5N3 <2642243996@qq.com>
Date: Sun, 16 Jan 2022 18:49:30 +0800
Subject: [PATCH 19/19] Inference fix

Update reducedim.jl
---
 base/reducedim.jl | 11 ++++++-----
 1 file changed, 6 insertions(+), 5 deletions(-)

diff --git a/base/reducedim.jl b/base/reducedim.jl
index 1b437c556d4f3..d55db2768e62b 100644
--- a/base/reducedim.jl
+++ b/base/reducedim.jl
@@ -181,22 +181,23 @@ function reducedim_init(f::ExtremaMap, op::typeof(_extrema_rf), A::AbstractArray
     v0 = reverse(mapreduce(f, op, A1)) # turn minmax to maxmin
 
     T = _realtype(f.f, promote_union(eltype(A)))
-    Tr = v0[1] isa T && v0[2] isa T ? NTuple{2,T} : typeof(v0)
+    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
-        T1, T2 = Tr.parameters
         # TODO: Some types, like BigInt, don't support typemin/typemax.
         # So a Matrix{Union{BigInt, Missing}} can still error here.
-        v0 = typemax(nonmissingtype(T1)), typemin(nonmissingtype(T2))
+        v0 = typemax(nonmissingtype(Tmin)), typemin(nonmissingtype(Tmax))
     end
     # v0 may have changed type.
-    Tr = v0[1] isa T && v0[2] isa T ? NTuple{2,T} : typeof(v0)
+    Tmin = v0[1] isa T ? T : typeof(v0[1])
+    Tmax = v0[2] isa T ? T : typeof(v0[2])
 
-    return reducedim_initarray(A, region, v0, Tr)
+    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} =