From 34f8126ffd11bf19b13ccd1a4b78668e1f831579 Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Thu, 14 Mar 2024 15:40:48 -0400 Subject: [PATCH 1/3] fix literal-pow to return the right type when the base is -1 --- base/intfuncs.jl | 19 ++++++++++++++++--- 1 file changed, 16 insertions(+), 3 deletions(-) diff --git a/base/intfuncs.jl b/base/intfuncs.jl index 916cd28a32c6f..be2f1a218a51e 100644 --- a/base/intfuncs.jl +++ b/base/intfuncs.jl @@ -369,8 +369,19 @@ const HWNumber = Union{HWReal, Complex{<:HWReal}, Rational{<:HWReal}} @inline literal_pow(::typeof(^), x::HWNumber, ::Val{1}) = x @inline literal_pow(::typeof(^), x::HWNumber, ::Val{2}) = x*x @inline literal_pow(::typeof(^), x::HWNumber, ::Val{3}) = x*x*x -@inline literal_pow(::typeof(^), x::HWNumber, ::Val{-1}) = inv(x) -@inline literal_pow(::typeof(^), x::HWNumber, ::Val{-2}) = (i=inv(x); i*i) +@inline function literal_pow(::typeof(^), x::HWNumber, ::Val{-1}) + if x == -one(x) + return isodd(p) ? one(x) : -one(x) + end + inv(x) +end +@inline function literal_pow(::typeof(^), x::HWNumber, ::Val{-2}) + if x == -one(x) + return isodd(p) ? one(x) : -one(x) + end + i=inv(x) + return i*i +end # don't use the inv(x) transformation here since float^p is slightly more accurate @inline literal_pow(::typeof(^), x::AbstractFloat, ::Val{p}) where {p} = x^p @@ -379,7 +390,9 @@ const HWNumber = Union{HWReal, Complex{<:HWReal}, Rational{<:HWReal}} # for other types, define x^-n as inv(x)^n so that negative literal powers can # be computed in a type-stable way even for e.g. integers. @inline function literal_pow(f::typeof(^), x, ::Val{p}) where {p} - if p < 0 + if x == -one(x) + return isodd(p) ? one(x) : -one(x) + elseif p < 0 if x isa BitInteger64 f(Float64(x), p) # inv would cause rounding, while Float64^Integer is able to compensate the inverse else From a5f8b2a0b7ed69a99dbd5adae9ab21e3d7a80ba5 Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Thu, 14 Mar 2024 15:44:03 -0400 Subject: [PATCH 2/3] add tests --- test/math.jl | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/test/math.jl b/test/math.jl index f551bb3a5d4b7..eeed7b5f0f916 100644 --- a/test/math.jl +++ b/test/math.jl @@ -1459,6 +1459,13 @@ end # two cases where we have observed > 1 ULP in the past @test 0.0013653274095082324^-97.60372292227069 == 4.088393948750035e279 @test 8.758520413376658e-5^70.55863059215994 == 5.052076767078296e-287 + # test literal pow for base = -1 + @test (-1)^(-1) === (-1)^((-1,)[1]) + @test (-1)^(-2) === (-1)^((-2,)[1]) + @test (-1)^(-3) === (-1)^((-3,)[1]) + @test (-1)^(1) === (-1)^((1,)[1]) + @test (-1)^(2) === (-1)^((2,)[1]) + @test (-1)^(3) === (-1)^((3,)[1]) end # Test that sqrt behaves correctly and doesn't exhibit fp80 double rounding. From db3677b505da297f6db6850aff2d241ce646969d Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Fri, 15 Mar 2024 09:50:30 -0400 Subject: [PATCH 3/3] Update base/intfuncs.jl Co-authored-by: Max Horn --- base/intfuncs.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/base/intfuncs.jl b/base/intfuncs.jl index be2f1a218a51e..7a3c9d3ed7816 100644 --- a/base/intfuncs.jl +++ b/base/intfuncs.jl @@ -371,7 +371,7 @@ const HWNumber = Union{HWReal, Complex{<:HWReal}, Rational{<:HWReal}} @inline literal_pow(::typeof(^), x::HWNumber, ::Val{3}) = x*x*x @inline function literal_pow(::typeof(^), x::HWNumber, ::Val{-1}) if x == -one(x) - return isodd(p) ? one(x) : -one(x) + return x end inv(x) end