From 1a4640cc4b25c128eeb09692223ca65c2a38dff2 Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Fri, 19 Aug 2022 11:38:35 -0400 Subject: [PATCH 1/4] fixes --- base/math.jl | 17 ++++++++++------- 1 file changed, 10 insertions(+), 7 deletions(-) diff --git a/base/math.jl b/base/math.jl index 8a4735c43ebea..64f9e9e1a78d4 100644 --- a/base/math.jl +++ b/base/math.jl @@ -1098,14 +1098,17 @@ end # @constprop aggressive to help the compiler see the switch between the integer and float # variants for callers with constant `y` @constprop :aggressive function ^(x::Float64, y::Float64) - yint = unsafe_trunc(Int, y) # Note, this is actually safe since julia freezes the result - y == yint && return x^yint - #numbers greater than 2*inv(eps(T)) must be even, and the pow will overflow - y >= 2*inv(eps()) && return x^(typemax(Int64)-1) + # Exponents greater than this will always overflow or underflow. + # Note that NaN can pass through this, but that will end up fine. + if abs(y)>0x1.8p62 + y = sign(y)*0x1.8p62 + end + yint = unsafe_trunc(Int, y) # This is actually safe since julia freezes the result + y == yint && return @noinline x^yint xu = reinterpret(UInt64, x) - x<0 && y > -4e18 && throw_exp_domainerror(x) # |y| is small enough that y isn't an integer - x === 1.0 && return 1.0 - x==0 && return abs(y)*Inf*(!(y>0)) + x<0 && throw_exp_domainerror(x) # |y| is small enough that y isn't an integer + xu == reinterpret(UInt64, 1.0) && return 1.0 + 2*xu==0 && return abs(y)*Inf*(!(y>0)) # if x==0 !isfinite(x) && return x*(y>0 || isnan(x)) # x is inf or NaN if xu < (UInt64(1)<<52) # x is subnormal xu = reinterpret(UInt64, x * 0x1p52) # normalize x From 21d2d24e908b680273476c8ddf62d4ce43cf064d Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Fri, 19 Aug 2022 12:36:40 -0400 Subject: [PATCH 2/4] fix tests and code --- base/math.jl | 28 +++++++++++++++++----------- test/math.jl | 8 ++++++-- 2 files changed, 23 insertions(+), 13 deletions(-) diff --git a/base/math.jl b/base/math.jl index 64f9e9e1a78d4..e2a135c339acb 100644 --- a/base/math.jl +++ b/base/math.jl @@ -1098,17 +1098,18 @@ end # @constprop aggressive to help the compiler see the switch between the integer and float # variants for callers with constant `y` @constprop :aggressive function ^(x::Float64, y::Float64) + xu = reinterpret(UInt64, x) + xu == reinterpret(UInt64, 1.0) && return 1.0 # Exponents greater than this will always overflow or underflow. # Note that NaN can pass through this, but that will end up fine. - if abs(y)>0x1.8p62 + if !(abs(y)<0x1.8p62) + isnan(y) && return y y = sign(y)*0x1.8p62 end yint = unsafe_trunc(Int, y) # This is actually safe since julia freezes the result y == yint && return @noinline x^yint - xu = reinterpret(UInt64, x) - x<0 && throw_exp_domainerror(x) # |y| is small enough that y isn't an integer - xu == reinterpret(UInt64, 1.0) && return 1.0 2*xu==0 && return abs(y)*Inf*(!(y>0)) # if x==0 + x<0 && throw_exp_domainerror(x) # |y| is small enough that y isn't an integer !isfinite(x) && return x*(y>0 || isnan(x)) # x is inf or NaN if xu < (UInt64(1)<<52) # x is subnormal xu = reinterpret(UInt64, x * 0x1p52) # normalize x @@ -1127,18 +1128,23 @@ end end @constprop :aggressive function ^(x::T, y::T) where T <: Union{Float16, Float32} - yint = unsafe_trunc(Int64, y) # Note, this is actually safe since julia freezes the result + x == 1 && return one(T) + # Exponents greater than this will always overflow or underflow. + # Note that NaN can pass through this, but that will end up fine. + max_exp = T == Float16 ? T(3<<14) : T(0x1.Ap30) + if !(abs(y)= 2*inv(eps(T)) && return x^(typemax(Int64)-1) - x < 0 && y > -4e18 && throw_exp_domainerror(x) # |y| is small enough that y isn't an integer + x < 0 && throw_exp_domainerror(x) # |y| is small enough that y isn't an integer + !isfinite(x) && return x*(y>0 || isnan(x)) + x==0 && return abs(y)*T(Inf)*(!(y>0)) return pow_body(x, y) end @inline function pow_body(x::T, y::T) where T <: Union{Float16, Float32} - x == 1 && return one(T) - !isfinite(x) && return x*(y>0 || isnan(x)) - x==0 && return abs(y)*T(Inf)*(!(y>0)) return T(exp2(log2(abs(widen(x))) * y)) end diff --git a/test/math.jl b/test/math.jl index bae1f571ef16a..6ba40b7daa968 100644 --- a/test/math.jl +++ b/test/math.jl @@ -1325,8 +1325,12 @@ end for T in (Float16, Float32, Float64) for x in (0.0, -0.0, 1.0, 10.0, 2.0, Inf, NaN, -Inf, -NaN) for y in (0.0, -0.0, 1.0, -3.0,-10.0 , Inf, NaN, -Inf, -NaN) - got, expected = T(x)^T(y), T(big(x))^T(y) - @test isnan_type(T, got) && isnan_type(T, expected) || (got === expected) + got, expected = T(x)^T(y), T(big(x)^T(y)) + if isnan(expected) + @test isnan_type(T, got) || T.((x,y)) + else + @test got == expected || T.((x,y)) + end end end for _ in 1:2^16 From e149aaaef71615a5418e0cb8482f2b606b53f73e Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Fri, 19 Aug 2022 14:31:54 -0400 Subject: [PATCH 3/4] Int!=Int64 --- base/math.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/base/math.jl b/base/math.jl index e2a135c339acb..d0c4ff0fef882 100644 --- a/base/math.jl +++ b/base/math.jl @@ -1106,7 +1106,7 @@ end isnan(y) && return y y = sign(y)*0x1.8p62 end - yint = unsafe_trunc(Int, y) # This is actually safe since julia freezes the result + yint = unsafe_trunc(Int64, y) # This is actually safe since julia freezes the result y == yint && return @noinline x^yint 2*xu==0 && return abs(y)*Inf*(!(y>0)) # if x==0 x<0 && throw_exp_domainerror(x) # |y| is small enough that y isn't an integer From 6857a3dd832d56eff8bdc799e027cb9a71616802 Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Fri, 19 Aug 2022 14:32:24 -0400 Subject: [PATCH 4/4] remove un-needed comment --- base/math.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/base/math.jl b/base/math.jl index d0c4ff0fef882..f1ee129305418 100644 --- a/base/math.jl +++ b/base/math.jl @@ -1138,7 +1138,7 @@ end end yint = unsafe_trunc(Int32, y) # This is actually safe since julia freezes the result y == yint && return x^yint - x < 0 && throw_exp_domainerror(x) # |y| is small enough that y isn't an integer + x < 0 && throw_exp_domainerror(x) !isfinite(x) && return x*(y>0 || isnan(x)) x==0 && return abs(y)*T(Inf)*(!(y>0)) return pow_body(x, y)