From 8ed3bc6fc4cb1654cd45654f60091c5f9521e7fd Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Mon, 22 Aug 2022 04:12:28 -0400 Subject: [PATCH] backport #46412 to 1.8 (#46417) --- base/math.jl | 37 +++++++++++++++++++++++-------------- test/math.jl | 29 ++++++++++++++++++++++++++--- 2 files changed, 49 insertions(+), 17 deletions(-) diff --git a/base/math.jl b/base/math.jl index 1f194f73ca7f9..eb69c555d985b 100644 --- a/base/math.jl +++ b/base/math.jl @@ -996,18 +996,22 @@ 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) - x<0 && y > -4e18 && throw_exp_domainerror(x) # |y| is small enough that y isn't an integer x == 1 && 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) + isnan(y) && return y + y = sign(y)*0x1.8p62 + end + yint = unsafe_trunc(Int64, y) # This is actually safe since julia freezes the result + y == yint && return x^yint + x<0 && throw_exp_domainerror(x) + !isfinite(x) && return x*(y>0 || isnan(x)) + x==0 && return abs(y)*Inf*(!(y>0)) return pow_body(x, y) end @inline function pow_body(x::Float64, y::Float64) - !isfinite(x) && return x*(y>0 || isnan(x)) - x==0 && return abs(y)*Inf*(!(y>0)) logxhi,logxlo = Base.Math._log_ext(x) xyhi, xylo = two_mul(logxhi,y) xylo = muladd(logxlo, y, xylo) @@ -1016,18 +1020,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) + !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 5daf1ef465509..56d97562ca0ab 100644 --- a/test/math.jl +++ b/test/math.jl @@ -1318,18 +1318,41 @@ end end @testset "pow" begin + # tolerance by type for regular powers + POW_TOLS = Dict(Float16=>[.51, .51, 2.0, 1.5], + Float32=>[.51, .51, 2.0, 1.5], + Float64=>[1.0, 1.5, 2.0, 1.5]) 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 + # note x won't be subnormal here x=rand(T)*100; y=rand(T)*200-100 got, expected = x^y, widen(x)^y if isfinite(eps(T(expected))) - @test abs(expected-got) <= 1.3*eps(T(expected)) || (x,y) + if y == T(-2) # unfortunately x^-2 is less accurate for performance reasons. + @test abs(expected-got) <= POW_TOLS[T][3]*eps(T(expected)) || (x,y) + elseif y == T(3) # unfortunately x^3 is less accurate for performance reasons. + @test abs(expected-got) <= POW_TOLS[T][4]*eps(T(expected)) || (x,y) + else + @test abs(expected-got) <= POW_TOLS[T][1]*eps(T(expected)) || (x,y) + end + end + end + for _ in 1:2^14 + # test subnormal(x), y in -1.2, 1.8 since anything larger just overflows. + x=rand(T)*floatmin(T); y=rand(T)*3-T(1.2) + got, expected = x^y, widen(x)^y + if isfinite(eps(T(expected))) + @test abs(expected-got) <= POW_TOLS[T][2]*eps(T(expected)) || (x,y) end end # test (-x)^y for y larger than typemax(Int)