From 59365745d8f669628d2e06c520ed4246e6f490d7 Mon Sep 17 00:00:00 2001 From: N5N3 <2642243996@qq.com> Date: Tue, 17 Aug 2021 23:43:31 +0800 Subject: [PATCH 1/6] Fast `min`/`max/minmax` for `Float32/64` https://discourse.julialang.org/t/faster-min-max-for-float64-32/45464/40 --- base/math.jl | 17 ++++++++--------- 1 file changed, 8 insertions(+), 9 deletions(-) diff --git a/base/math.jl b/base/math.jl index 9550a0a54b496..9cc44f6c9c4f9 100644 --- a/base/math.jl +++ b/base/math.jl @@ -758,17 +758,16 @@ end atan(y::Real, x::Real) = atan(promote(float(y),float(x))...) atan(y::T, x::T) where {T<:AbstractFloat} = Base.no_op_err("atan", T) -max(x::T, y::T) where {T<:AbstractFloat} = ifelse((y > x) | (signbit(y) < signbit(x)), - ifelse(isnan(x), x, y), ifelse(isnan(y), y, x)) +_isless(x::T, y::T) where {T<:AbstractFloat} = (x < y) || (signbit(x) > signbit(y)) +min(x::T, y::T) where {T<:AbstractFloat} = isnan(x) || ~isnan(y) && _isless(x, y) ? x : y +max(x::T, y::T) where {T<:AbstractFloat} = isnan(x) || ~isnan(y) && _isless(y, x) ? x : y +minmax(x::T, y::T) where {T<:AbstractFloat} = min(x, y), max(x, y) +_isless(x::T, y::T) where {T<:Union{Float32,Float64}} = signbit(x - y) +min(x::T, y::T) where {T<:Union{Float32,Float64}} = ifelse(isnan(x) | ~isnan(y) & _isless(x, y), x, y) +max(x::T, y::T) where {T<:Union{Float32,Float64}} = ifelse(isnan(x) | ~isnan(y) & _isless(y, x), x, y) -min(x::T, y::T) where {T<:AbstractFloat} = ifelse((y < x) | (signbit(y) > signbit(x)), - ifelse(isnan(x), x, y), ifelse(isnan(y), y, x)) - -minmax(x::T, y::T) where {T<:AbstractFloat} = - ifelse(isnan(x) | isnan(y), ifelse(isnan(x), (x,x), (y,y)), - ifelse((y > x) | (signbit(x) > signbit(y)), (x,y), (y,x))) - +_isless(x::Float16, y::Float16) = _isless(widen(x), widen(y)) """ ldexp(x, n) From 3e39c57a051da7b46897cee1b8c73b800b613477 Mon Sep 17 00:00:00 2001 From: N5N3 <2642243996@qq.com> Date: Fri, 21 Jan 2022 09:33:22 +0800 Subject: [PATCH 2/6] extrema optimization for IEEEFloat --- base/reduce.jl | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/base/reduce.jl b/base/reduce.jl index 45284d884a279..2b987b7b9dff6 100644 --- a/base/reduce.jl +++ b/base/reduce.jl @@ -847,8 +847,15 @@ 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 IEEEFloat +function _extrema_rf(x::NTuple{2,T}, y::NTuple{2,T}) where {T<:IEEEFloat} + (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 From 04f07fb0a130e491d8c71fe7ba0ac03fe0ea11fc Mon Sep 17 00:00:00 2001 From: N5N3 <2642243996@qq.com> Date: Fri, 21 Jan 2022 10:37:43 +0800 Subject: [PATCH 3/6] Add BigFloat related optimization Update mpfr.jl --- base/mpfr.jl | 27 ++++++++++++++------------- 1 file changed, 14 insertions(+), 13 deletions(-) diff --git a/base/mpfr.jl b/base/mpfr.jl index 60f59cdb0af7e..97e4535c065d2 100644 --- a/base/mpfr.jl +++ b/base/mpfr.jl @@ -16,7 +16,7 @@ import cosh, sinh, tanh, sech, csch, coth, acosh, asinh, atanh, lerpi, cbrt, typemax, typemin, unsafe_trunc, floatmin, floatmax, rounding, setrounding, maxintfloat, widen, significand, frexp, tryparse, iszero, - isone, big, _string_n, decompose + isone, big, _string_n, decompose, minmax import ..Rounding: rounding_raw, setrounding_raw @@ -697,20 +697,21 @@ function log1p(x::BigFloat) return z end -function max(x::BigFloat, y::BigFloat) - isnan(x) && return x - isnan(y) && return y - z = BigFloat() - ccall((:mpfr_max, :libmpfr), Int32, (Ref{BigFloat}, Ref{BigFloat}, Ref{BigFloat}, MPFRRoundingMode), z, x, y, ROUNDING_MODE[]) - return z +# For `min`/`max`, general fallback for `AbstractFloat` is good enough. +# Only implement `minmax` and `_extrema_rf` to avoid repeated calls. +function minmax(x::BigFloat, y::BigFloat) + isnan(x) && return x, x + isnan(y) && return y, y + Base.Math._isless(x, y) ? (x, y) : (y, x) end -function min(x::BigFloat, y::BigFloat) - isnan(x) && return x - isnan(y) && return y - z = BigFloat() - ccall((:mpfr_min, :libmpfr), Int32, (Ref{BigFloat}, Ref{BigFloat}, Ref{BigFloat}, MPFRRoundingMode), z, x, y, ROUNDING_MODE[]) - return z +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 = Base.Math._isless(x1, y1) ? x1 : y1 + z2 = Base.Math._isless(x2, y2) ? y2 : x2 + z1, z2 end function modf(x::BigFloat) From 39d68783cdf142ce333fa22596da04437128cc1d Mon Sep 17 00:00:00 2001 From: N5N3 <2642243996@qq.com> Date: Fri, 21 Jan 2022 10:37:11 +0800 Subject: [PATCH 4/6] Test `min`/`max`/`minmax` for all float types --- test/numbers.jl | 66 ++++++++++++++++++++++++++++--------------------- 1 file changed, 38 insertions(+), 28 deletions(-) diff --git a/test/numbers.jl b/test/numbers.jl index ad521d7382713..2d55f37de643e 100644 --- a/test/numbers.jl +++ b/test/numbers.jl @@ -95,34 +95,44 @@ end @test max(1) === 1 @test minmax(1) === (1, 1) @test minmax(5, 3) == (3, 5) - @test minmax(3., 5.) == (3., 5.) - @test minmax(5., 3.) == (3., 5.) - @test minmax(3., NaN) ≣ (NaN, NaN) - @test minmax(NaN, 3) ≣ (NaN, NaN) - @test minmax(Inf, NaN) ≣ (NaN, NaN) - @test minmax(NaN, Inf) ≣ (NaN, NaN) - @test minmax(-Inf, NaN) ≣ (NaN, NaN) - @test minmax(NaN, -Inf) ≣ (NaN, NaN) - @test minmax(NaN, NaN) ≣ (NaN, NaN) - @test min(-0.0,0.0) === min(0.0,-0.0) - @test max(-0.0,0.0) === max(0.0,-0.0) - @test minmax(-0.0,0.0) === minmax(0.0,-0.0) - @test max(-3.2, 5.1) == max(5.1, -3.2) == 5.1 - @test min(-3.2, 5.1) == min(5.1, -3.2) == -3.2 - @test max(-3.2, Inf) == max(Inf, -3.2) == Inf - @test max(-3.2, NaN) ≣ max(NaN, -3.2) ≣ NaN - @test min(5.1, Inf) == min(Inf, 5.1) == 5.1 - @test min(5.1, -Inf) == min(-Inf, 5.1) == -Inf - @test min(5.1, NaN) ≣ min(NaN, 5.1) ≣ NaN - @test min(5.1, -NaN) ≣ min(-NaN, 5.1) ≣ NaN - @test minmax(-3.2, 5.1) == (min(-3.2, 5.1), max(-3.2, 5.1)) - @test minmax(-3.2, Inf) == (min(-3.2, Inf), max(-3.2, Inf)) - @test minmax(-3.2, NaN) ≣ (min(-3.2, NaN), max(-3.2, NaN)) - @test (max(Inf,NaN), max(-Inf,NaN), max(Inf,-NaN), max(-Inf,-NaN)) ≣ (NaN,NaN,NaN,NaN) - @test (max(NaN,Inf), max(NaN,-Inf), max(-NaN,Inf), max(-NaN,-Inf)) ≣ (NaN,NaN,NaN,NaN) - @test (min(Inf,NaN), min(-Inf,NaN), min(Inf,-NaN), min(-Inf,-NaN)) ≣ (NaN,NaN,NaN,NaN) - @test (min(NaN,Inf), min(NaN,-Inf), min(-NaN,Inf), min(-NaN,-Inf)) ≣ (NaN,NaN,NaN,NaN) - @test minmax(-Inf,NaN) ≣ (min(-Inf,NaN), max(-Inf,NaN)) + Top(T, op, x, y) = op(T.(x), T.(y)) + Top(T, op) = (x, y) -> Top(T, op, x, y) + _compare(x, y) = x == y + for T in (Float16, Float32, Float64, BigFloat) + minmax = Top(T,Base.minmax) + min = Top(T,Base.min) + max = Top(T,Base.max) + (==) = Top(T,_compare) + (===) = Top(T,Base.isequal) # we only use === to compare -0.0/0.0, `isequal` should be equalvient + @test minmax(3., 5.) == (3., 5.) + @test minmax(5., 3.) == (3., 5.) + @test minmax(3., NaN) ≣ (NaN, NaN) + @test minmax(NaN, 3) ≣ (NaN, NaN) + @test minmax(Inf, NaN) ≣ (NaN, NaN) + @test minmax(NaN, Inf) ≣ (NaN, NaN) + @test minmax(-Inf, NaN) ≣ (NaN, NaN) + @test minmax(NaN, -Inf) ≣ (NaN, NaN) + @test minmax(NaN, NaN) ≣ (NaN, NaN) + @test min(-0.0,0.0) === min(0.0,-0.0) + @test max(-0.0,0.0) === max(0.0,-0.0) + @test minmax(-0.0,0.0) === minmax(0.0,-0.0) + @test max(-3.2, 5.1) == max(5.1, -3.2) == 5.1 + @test min(-3.2, 5.1) == min(5.1, -3.2) == -3.2 + @test max(-3.2, Inf) == max(Inf, -3.2) == Inf + @test max(-3.2, NaN) ≣ max(NaN, -3.2) ≣ NaN + @test min(5.1, Inf) == min(Inf, 5.1) == 5.1 + @test min(5.1, -Inf) == min(-Inf, 5.1) == -Inf + @test min(5.1, NaN) ≣ min(NaN, 5.1) ≣ NaN + @test min(5.1, -NaN) ≣ min(-NaN, 5.1) ≣ NaN + @test minmax(-3.2, 5.1) == (min(-3.2, 5.1), max(-3.2, 5.1)) + @test minmax(-3.2, Inf) == (min(-3.2, Inf), max(-3.2, Inf)) + @test minmax(-3.2, NaN) ≣ (min(-3.2, NaN), max(-3.2, NaN)) + @test (max(Inf,NaN), max(-Inf,NaN), max(Inf,-NaN), max(-Inf,-NaN)) ≣ (NaN,NaN,NaN,NaN) + @test (max(NaN,Inf), max(NaN,-Inf), max(-NaN,Inf), max(-NaN,-Inf)) ≣ (NaN,NaN,NaN,NaN) + @test (min(Inf,NaN), min(-Inf,NaN), min(Inf,-NaN), min(-Inf,-NaN)) ≣ (NaN,NaN,NaN,NaN) + @test (min(NaN,Inf), min(NaN,-Inf), min(-NaN,Inf), min(-NaN,-Inf)) ≣ (NaN,NaN,NaN,NaN) + @test minmax(-Inf,NaN) ≣ (min(-Inf,NaN), max(-Inf,NaN)) + end end @testset "fma" begin let x = Int64(7)^7 From 18f5a7c4d3171b2f147b45465177c2ae0f1e74f2 Mon Sep 17 00:00:00 2001 From: N5N3 <2642243996@qq.com> Date: Fri, 21 Jan 2022 10:55:38 +0800 Subject: [PATCH 5/6] add _extrema_rf test --- test/numbers.jl | 24 ++++++++++++++++++++++++ 1 file changed, 24 insertions(+) diff --git a/test/numbers.jl b/test/numbers.jl index 2d55f37de643e..e972064fb8f10 100644 --- a/test/numbers.jl +++ b/test/numbers.jl @@ -134,6 +134,30 @@ end @test minmax(-Inf,NaN) ≣ (min(-Inf,NaN), max(-Inf,NaN)) end end +@testset "Base._extrema_rf for float" begin + for T in (Float16, Float32, Float64, BigFloat) + ordered = T[-Inf, -5, -0.0, 0.0, 3, Inf] + unorded = T[NaN, -NaN] + for i1 in 1:6, i2 in 1:6, j1 in 1:6, j2 in 1:6 + x = ordered[i1], ordered[i2] + y = ordered[j1], ordered[j2] + z = ordered[min(i1,j1)], ordered[max(i2,j2)] + @test Base._extrema_rf(x, y) === z + end + for i in 1:2, j1 in 1:6, j2 in 1:6 # unordered test (only 1 NaN) + x = unorded[i] , unorded[i] + y = ordered[j1], ordered[j2] + @test Base._extrema_rf(x, y) === x + @test Base._extrema_rf(y, x) === x + end + for i in 1:2, j in 1:2 # unordered test (2 NaNs) + x = unorded[i], unorded[i] + y = unorded[j], unorded[j] + z = Base._extrema_rf(x, y) + @test z === x || z === y + end + end +end @testset "fma" begin let x = Int64(7)^7 @test fma(x-1, x-2, x-3) == (x-1) * (x-2) + (x-3) From b64fddb2db42d41af1f3e1ca61e1785406183676 Mon Sep 17 00:00:00 2001 From: N5N3 <2642243996@qq.com> Date: Wed, 1 Jun 2022 09:58:50 +0800 Subject: [PATCH 6/6] Some performance tune. --- base/math.jl | 26 ++++++++++++++++++++++---- 1 file changed, 22 insertions(+), 4 deletions(-) diff --git a/base/math.jl b/base/math.jl index 9cc44f6c9c4f9..81a3c453ad88e 100644 --- a/base/math.jl +++ b/base/math.jl @@ -763,11 +763,29 @@ min(x::T, y::T) where {T<:AbstractFloat} = isnan(x) || ~isnan(y) && _isless(x, y max(x::T, y::T) where {T<:AbstractFloat} = isnan(x) || ~isnan(y) && _isless(y, x) ? x : y minmax(x::T, y::T) where {T<:AbstractFloat} = min(x, y), max(x, y) -_isless(x::T, y::T) where {T<:Union{Float32,Float64}} = signbit(x - y) -min(x::T, y::T) where {T<:Union{Float32,Float64}} = ifelse(isnan(x) | ~isnan(y) & _isless(x, y), x, y) -max(x::T, y::T) where {T<:Union{Float32,Float64}} = ifelse(isnan(x) | ~isnan(y) & _isless(y, x), x, y) +_isless(x::Float16, y::Float16) = signbit(widen(x) - widen(y)) -_isless(x::Float16, y::Float16) = _isless(widen(x), widen(y)) +function min(x::T, y::T) where {T<:Union{Float32,Float64}} + diff = x - y + argmin = ifelse(signbit(diff), x, y) + anynan = isnan(x)|isnan(y) + ifelse(anynan, diff, argmin) +end + +function max(x::T, y::T) where {T<:Union{Float32,Float64}} + diff = x - y + argmax = ifelse(signbit(diff), y, x) + anynan = isnan(x)|isnan(y) + ifelse(anynan, diff, argmax) +end + +function minmax(x::T, y::T) where {T<:Union{Float32,Float64}} + diff = x - y + sdiff = signbit(diff) + min, max = ifelse(sdiff, x, y), ifelse(sdiff, y, x) + anynan = isnan(x)|isnan(y) + ifelse(anynan, diff, min), ifelse(anynan, diff, max) +end """ ldexp(x, n)