improvements to accuracy/performance for float^integer (#24500)

JuliaLang · Dec 5, 2017 · 8fa23ed · 8fa23ed
1 parent 9fe9d2c
commit 8fa23ed
Show file tree

Hide file tree

Showing 3 changed files with 18 additions and 3 deletions.
diff --git a/base/intfuncs.jl b/base/intfuncs.jl
@@ -234,6 +234,12 @@ const HWNumber = Union{HWReal, Complex{<:HWReal}, Rational{<:HWReal}}
 @inline literal_pow(::typeof(^), x::HWNumber, ::Val{2}) = x*x
 @inline literal_pow(::typeof(^), x::HWNumber, ::Val{3}) = x*x*x
 
+# 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
+@inline literal_pow(::typeof(^), x::AbstractFloat, ::Val{-1}) = inv(x)
+
+# 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 @generated function literal_pow(f::typeof(^), x, ::Val{p}) where {p}
     if p < 0
         :(literal_pow(^, inv(x), $(Val{-p}())))

diff --git a/base/math.jl b/base/math.jl
@@ -732,9 +732,9 @@ end
     end
     z
 end
-@inline ^(x::Float64, y::Integer) = x ^ Float64(y)
-@inline ^(x::Float32, y::Integer) = x ^ Float32(y)
-@inline ^(x::Float16, y::Integer) = Float16(Float32(x) ^ Float32(y))
+@inline ^(x::Float64, y::Integer) = ccall("llvm.pow.f64", llvmcall, Float64, (Float64, Float64), x, Float64(y))
+@inline ^(x::Float32, y::Integer) = ccall("llvm.pow.f32", llvmcall, Float32, (Float32, Float32), x, Float32(y))
+@inline ^(x::Float16, y::Integer) = Float16(Float32(x) ^ y)
 @inline literal_pow(::typeof(^), x::Float16, ::Val{p}) where {p} = Float16(literal_pow(^,Float32(x),Val(p)))
 
 function angle_restrict_symm(theta)

diff --git a/test/numbers.jl b/test/numbers.jl
@@ -2965,6 +2965,15 @@ module M20889 # do we get the expected behavior without importing Base.^?
     Test.@test PR20889(2)^3 == 5
 end
 
+@testset "literal negative power accuracy" begin
+    # a few cases chosen to maximize the error for inv(x)^+n:
+    @test 0.7130409001548401^-2 === 1.9668494399322154
+    @test 0.09496527f0^-2 === 110.88438f0
+    @test 0.20675883960662367^-100 === 2.841786121808917e68
+    @test 0.6123676f0^-100 === 1.9896729f21
+    @test 0.004155780785470562^-1 === 240.6286692253353
+end
+
 @testset "iszero & isone" begin
     # Numeric scalars
     for T in [Float16, Float32, Float64, BigFloat,