From 14a46775b1b3d475a7619671d33fedad9ff78aac Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Mon, 6 Nov 2017 17:06:32 -0500 Subject: [PATCH 1/3] improvements to accuracy/performance for float^integer --- base/intfuncs.jl | 7 +++++++ base/math.jl | 6 +++--- base/mpfr.jl | 2 ++ test/math.jl | 8 ++++++++ 4 files changed, 20 insertions(+), 3 deletions(-) diff --git a/base/intfuncs.jl b/base/intfuncs.jl index cffb2323a13dc..b178447999080 100644 --- a/base/intfuncs.jl +++ b/base/intfuncs.jl @@ -234,6 +234,13 @@ 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 x^p is slightly more accurate +@inline literal_pow(::typeof(^), x::Union{Float32,Float64}, ::Val{p}) where {p} = x^p +@inline literal_pow(::typeof(^), x::Union{Float32,Float64}, ::Val{0}) = one(x) +@inline literal_pow(::typeof(^), x::Union{Float32,Float64}, ::Val{1}) = x +@inline literal_pow(::typeof(^), x::Union{Float32,Float64}, ::Val{2}) = x*x +@inline literal_pow(::typeof(^), x::Union{Float32,Float64}, ::Val{3}) = x*x*x + @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 index 9b4b317cc34ee..2c02fe5625c07 100644 --- 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/base/mpfr.jl b/base/mpfr.jl index 743a814e1af88..3c02ff8337f91 100644 --- a/base/mpfr.jl +++ b/base/mpfr.jl @@ -522,6 +522,8 @@ end ^(x::BigFloat, y::Integer) = typemin(Clong) <= y <= typemax(Clong) ? x^Clong(y) : x^BigInt(y) ^(x::BigFloat, y::Unsigned) = typemin(Culong) <= y <= typemax(Culong) ? x^Culong(y) : x^BigInt(y) +literal_pow(::typeof(^), x::BigFloat, ::Val{p}) where {p} = x^p + for f in (:exp, :exp2, :exp10, :expm1, :cosh, :sinh, :tanh, :sech, :csch, :coth, :cbrt) @eval function $f(x::BigFloat) z = BigFloat() diff --git a/test/math.jl b/test/math.jl index 43a4327e90e6c..20633b2b7dc62 100644 --- a/test/math.jl +++ b/test/math.jl @@ -830,3 +830,11 @@ end @test isnan_type(T, acos(T(NaN))) end 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 +end From 0db74f2660d9ccf5e69ac1828b8d36b6ab13e4c8 Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Mon, 6 Nov 2017 17:21:57 -0500 Subject: [PATCH 2/3] eliminate inv(x)^n fallback for all AbstractFloat types --- base/intfuncs.jl | 4 +++- base/mpfr.jl | 2 -- test/math.jl | 8 -------- test/numbers.jl | 8 ++++++++ 4 files changed, 11 insertions(+), 11 deletions(-) diff --git a/base/intfuncs.jl b/base/intfuncs.jl index b178447999080..7f364f9e810a4 100644 --- a/base/intfuncs.jl +++ b/base/intfuncs.jl @@ -235,7 +235,9 @@ const HWNumber = Union{HWReal, Complex{<:HWReal}, Rational{<:HWReal}} @inline literal_pow(::typeof(^), x::HWNumber, ::Val{3}) = x*x*x # don't use the inv(x) transformation here since x^p is slightly more accurate -@inline literal_pow(::typeof(^), x::Union{Float32,Float64}, ::Val{p}) where {p} = x^p +@inline literal_pow(::typeof(^), x::AbstractFloat, ::Val{p}) where {p} = x^p + +# eliminate method ambiguities: @inline literal_pow(::typeof(^), x::Union{Float32,Float64}, ::Val{0}) = one(x) @inline literal_pow(::typeof(^), x::Union{Float32,Float64}, ::Val{1}) = x @inline literal_pow(::typeof(^), x::Union{Float32,Float64}, ::Val{2}) = x*x diff --git a/base/mpfr.jl b/base/mpfr.jl index 3c02ff8337f91..743a814e1af88 100644 --- a/base/mpfr.jl +++ b/base/mpfr.jl @@ -522,8 +522,6 @@ end ^(x::BigFloat, y::Integer) = typemin(Clong) <= y <= typemax(Clong) ? x^Clong(y) : x^BigInt(y) ^(x::BigFloat, y::Unsigned) = typemin(Culong) <= y <= typemax(Culong) ? x^Culong(y) : x^BigInt(y) -literal_pow(::typeof(^), x::BigFloat, ::Val{p}) where {p} = x^p - for f in (:exp, :exp2, :exp10, :expm1, :cosh, :sinh, :tanh, :sech, :csch, :coth, :cbrt) @eval function $f(x::BigFloat) z = BigFloat() diff --git a/test/math.jl b/test/math.jl index 20633b2b7dc62..43a4327e90e6c 100644 --- a/test/math.jl +++ b/test/math.jl @@ -830,11 +830,3 @@ end @test isnan_type(T, acos(T(NaN))) end 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 -end diff --git a/test/numbers.jl b/test/numbers.jl index b610076b17fa5..a0370c66f4e1c 100644 --- a/test/numbers.jl +++ b/test/numbers.jl @@ -2965,6 +2965,14 @@ 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 +end + @testset "iszero & isone" begin # Numeric scalars for T in [Float16, Float32, Float64, BigFloat, From a9354df17afe600c10c68c61654c7bee6b5d7058 Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Tue, 14 Nov 2017 14:34:05 -0500 Subject: [PATCH 3/3] method ambiguity code no longer needed, add back float^-1 = inv(x) --- base/intfuncs.jl | 11 ++++------- test/numbers.jl | 1 + 2 files changed, 5 insertions(+), 7 deletions(-) diff --git a/base/intfuncs.jl b/base/intfuncs.jl index 7f364f9e810a4..8e7ab8543b6e6 100644 --- a/base/intfuncs.jl +++ b/base/intfuncs.jl @@ -234,15 +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 x^p is slightly more accurate +# 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) -# eliminate method ambiguities: -@inline literal_pow(::typeof(^), x::Union{Float32,Float64}, ::Val{0}) = one(x) -@inline literal_pow(::typeof(^), x::Union{Float32,Float64}, ::Val{1}) = x -@inline literal_pow(::typeof(^), x::Union{Float32,Float64}, ::Val{2}) = x*x -@inline literal_pow(::typeof(^), x::Union{Float32,Float64}, ::Val{3}) = x*x*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/test/numbers.jl b/test/numbers.jl index a0370c66f4e1c..395e2fdbf00c8 100644 --- a/test/numbers.jl +++ b/test/numbers.jl @@ -2971,6 +2971,7 @@ end @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