From 07ded3c4132ccf21d1123d4aa9192c38f64d3718 Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Fri, 5 Nov 2021 19:31:49 -0400 Subject: [PATCH 01/13] don't inline float64^float64 --- base/math.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/base/math.jl b/base/math.jl index c0a0cbafceec4..98f458b8acad4 100644 --- a/base/math.jl +++ b/base/math.jl @@ -939,7 +939,7 @@ function modf(x::T) where T<:IEEEFloat return (rx, ix) end -@inline function ^(x::Float64, y::Float64) +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 x<0 && y > -4e18 && throw_exp_domainerror(x) # |y| is small enough that y isn't an integer From 10c4a6bccca6cc6b47c98fb0071a7d755ab0f53a Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Mon, 8 Nov 2021 10:02:30 -0500 Subject: [PATCH 02/13] fix literal_pow to be more accurate for large negative bases --- base/intfuncs.jl | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/base/intfuncs.jl b/base/intfuncs.jl index ec57f7f80d809..0a60362a0ab14 100644 --- a/base/intfuncs.jl +++ b/base/intfuncs.jl @@ -319,6 +319,9 @@ 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 literal_pow(::typeof(^), x::HWNumber, ::Val{-3}) = i=inv(x); i*i*i # 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 @@ -326,9 +329,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} +@inline function literal_pow(f::typeof(^), x, ::Val{p}) where {p<:Integer} if p < 0 - literal_pow(^, inv(x), Val(-p)) + f(Float64(x),p) else f(x, p) end From bcf5e27ad4cb2d478f2f7f7c360952f929ca3286 Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Mon, 8 Nov 2021 10:03:08 -0500 Subject: [PATCH 03/13] fixup! fix literal_pow to be more accurate for large negative bases --- base/intfuncs.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/base/intfuncs.jl b/base/intfuncs.jl index 0a60362a0ab14..58944542be2ef 100644 --- a/base/intfuncs.jl +++ b/base/intfuncs.jl @@ -331,7 +331,7 @@ const HWNumber = Union{HWReal, Complex{<:HWReal}, Rational{<:HWReal}} # be computed in a type-stable way even for e.g. integers. @inline function literal_pow(f::typeof(^), x, ::Val{p}) where {p<:Integer} if p < 0 - f(Float64(x),p) + f(Float64(x), p) else f(x, p) end From 9a3c9db28984f35cb9822b9dcef3dd028086fbbb Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Mon, 8 Nov 2021 11:05:41 -0500 Subject: [PATCH 04/13] fix matrix^int --- base/intfuncs.jl | 13 ++++++++++--- 1 file changed, 10 insertions(+), 3 deletions(-) diff --git a/base/intfuncs.jl b/base/intfuncs.jl index 58944542be2ef..cca60b5210279 100644 --- a/base/intfuncs.jl +++ b/base/intfuncs.jl @@ -320,8 +320,8 @@ 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 @inline literal_pow(::typeof(^), x::HWNumber, ::Val{-1}) = inv(x) -@inline literal_pow(::typeof(^), x::HWNumber, ::Val{-2}) = i=inv(x); i*i -@inline literal_pow(::typeof(^), x::HWNumber, ::Val{-3}) = i=inv(x); i*i*i +@inline literal_pow(::typeof(^), x::HWNumber, ::Val{-2}) = (i=inv(x); i*i) +@inline literal_pow(::typeof(^), x::HWNumber, ::Val{-3}) = (i=inv(x); i*i*i) # 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 @@ -329,13 +329,20 @@ 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<:Integer} +@inline function literal_pow(f::typeof(^), x::Integer, ::Val{p}) where {p<:Integer} if p < 0 f(Float64(x), p) else f(x, p) end end +@inline function literal_pow(f::typeof(^), x, ::Val{p}) where {p<:Integer} + if p < 0 + f(inv(x), p) + else + f(x, p) + end +end # note: it is tempting to add optimized literal_pow(::typeof(^), x, ::Val{n}) # methods here for various n, but this easily leads to method ambiguities From 0e7e1602ff8943ce44e5e876beb5cb14047749e5 Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Mon, 8 Nov 2021 11:14:32 -0500 Subject: [PATCH 05/13] restrict special case to small ints --- base/intfuncs.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/base/intfuncs.jl b/base/intfuncs.jl index cca60b5210279..0c90f0d22d2c6 100644 --- a/base/intfuncs.jl +++ b/base/intfuncs.jl @@ -329,7 +329,7 @@ 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::Integer, ::Val{p}) where {p<:Integer} +@inline function literal_pow(f::typeof(^), x::BitInteger64, ::Val{p}) where {p<:Integer} if p < 0 f(Float64(x), p) else From e34a79d8c67db96f757263ebd4ed89ecc50e1d92 Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Mon, 8 Nov 2021 11:36:01 -0500 Subject: [PATCH 06/13] fix broken --- base/intfuncs.jl | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/base/intfuncs.jl b/base/intfuncs.jl index 0c90f0d22d2c6..b808678cc3cee 100644 --- a/base/intfuncs.jl +++ b/base/intfuncs.jl @@ -329,16 +329,16 @@ 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::BitInteger64, ::Val{p}) where {p<:Integer} +@inline function literal_pow(f::typeof(^), x::BitInteger64, ::Val{p}) where {p} if p < 0 f(Float64(x), p) else f(x, p) end end -@inline function literal_pow(f::typeof(^), x, ::Val{p}) where {p<:Integer} +@inline function literal_pow(f::typeof(^), x, ::Val{p}) where {p} if p < 0 - f(inv(x), p) + f(inv(x), -p) else f(x, p) end From e3704bb0062f9cd22913aa633aabcc77254b4de5 Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Mon, 8 Nov 2021 11:48:55 -0500 Subject: [PATCH 07/13] fix part 2 --- base/intfuncs.jl | 13 +++++-------- 1 file changed, 5 insertions(+), 8 deletions(-) diff --git a/base/intfuncs.jl b/base/intfuncs.jl index b808678cc3cee..beacb52444ad6 100644 --- a/base/intfuncs.jl +++ b/base/intfuncs.jl @@ -329,16 +329,13 @@ 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::BitInteger64, ::Val{p}) where {p} - if p < 0 - f(Float64(x), p) - else - f(x, p) - end -end @inline function literal_pow(f::typeof(^), x, ::Val{p}) where {p} if p < 0 - f(inv(x), -p) + if p isa BitInteger64 + f(Float64(x), p) + else + f(inv(x), -p) + end else f(x, p) end From 1ef574f9eb303c2d3d2d70bf33f525395ed41158 Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Mon, 8 Nov 2021 13:01:20 -0500 Subject: [PATCH 08/13] Update base/intfuncs.jl Co-authored-by: Kristoffer Carlsson --- base/intfuncs.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/base/intfuncs.jl b/base/intfuncs.jl index beacb52444ad6..1bec3758b2334 100644 --- a/base/intfuncs.jl +++ b/base/intfuncs.jl @@ -331,7 +331,7 @@ const HWNumber = Union{HWReal, Complex{<:HWReal}, Rational{<:HWReal}} # 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 p isa BitInteger64 + if x isa BitInteger64 f(Float64(x), p) else f(inv(x), -p) From be196d7c101fd1ce2e921148e66be4116f517a3b Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Mon, 8 Nov 2021 13:19:54 -0500 Subject: [PATCH 09/13] why is this so annoying --- base/math.jl | 2 +- base/mathconstants.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/base/math.jl b/base/math.jl index 98f458b8acad4..a07fbddf6b293 100644 --- a/base/math.jl +++ b/base/math.jl @@ -963,7 +963,7 @@ end end # compensated power by squaring -@inline function ^(x::Float64, n::Integer) +function ^(x::Float64, n::Integer) n == 0 && return one(x) y = 1.0 xnlo = ynlo = 0.0 diff --git a/base/mathconstants.jl b/base/mathconstants.jl index 7e0dbda3c19a4..3bb4bb52ad07f 100644 --- a/base/mathconstants.jl +++ b/base/mathconstants.jl @@ -76,7 +76,7 @@ julia> Base.MathConstants.eulergamma julia> dx = 10^-6; julia> sum(-exp(-x) * log(x) for x in dx:dx:100) * dx -0.5772078382499134 +0.5772078382499133 ``` """ γ, const eulergamma = γ From ac3677f73b8ce94c272c83404a11a6240923a37f Mon Sep 17 00:00:00 2001 From: Kristoffer Carlsson Date: Mon, 8 Nov 2021 19:37:54 +0100 Subject: [PATCH 10/13] remove the last inlines. --- base/math.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/base/math.jl b/base/math.jl index a07fbddf6b293..228d24c031579 100644 --- a/base/math.jl +++ b/base/math.jl @@ -952,7 +952,7 @@ function ^(x::Float64, y::Float64) hi = xyhi+xylo return Base.Math.exp_impl(hi, xylo-(hi-xyhi), Val(:ℯ)) end -@inline function ^(x::T, y::T) where T <: Union{Float16, Float32} +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 y == yint && return x^yint x < 0 && y > -4e18 && throw_exp_domainerror(x) # |y| is small enough that y isn't an integer @@ -973,7 +973,7 @@ function ^(x::Float64, n::Integer) x = rx n = -n end - n==3 && return x*x*x #keep compatability with literal_pow + n == 3 && return x*x*x # keep compatibility with literal_pow while n > 1 if n&1 > 0 yn = x*y @@ -988,9 +988,9 @@ function ^(x::Float64, n::Integer) !isfinite(x) && return x*y return muladd(x, y, muladd(y, xnlo, x*ynlo)) end -@inline function ^(x::Float32, n::Integer) +function ^(x::Float32, n::Integer) n < 0 && return inv(x)^(-n) - n==3 && return x*x*x #keep compatability with literal_pow + n == 3 && return x*x*x #keep compatibility with literal_pow Float32(Base.power_by_squaring(Float64(x),n)) end @inline ^(x::Float16, y::Integer) = Float16(Float32(x) ^ y) From ce5b5f33f0de8023c9e65ad6dd0c30c9ef3b89a9 Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Mon, 8 Nov 2021 14:54:30 -0500 Subject: [PATCH 11/13] fix another bug --- base/math.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/base/math.jl b/base/math.jl index 228d24c031579..e549a519356b0 100644 --- a/base/math.jl +++ b/base/math.jl @@ -969,6 +969,7 @@ function ^(x::Float64, n::Integer) xnlo = ynlo = 0.0 if n < 0 rx = inv(x) + n==-2 && return x*x #keep compatability with literal_pow isfinite(x) && (xnlo = -fma(x, rx, -1.) * rx) x = rx n = -n From 0d5a68c610f3bdac427bd12d22082fe479f51322 Mon Sep 17 00:00:00 2001 From: oscarddssmith Date: Mon, 8 Nov 2021 16:18:07 -0500 Subject: [PATCH 12/13] inv(x)!=x --- base/math.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/base/math.jl b/base/math.jl index e549a519356b0..fe6cc1839a988 100644 --- a/base/math.jl +++ b/base/math.jl @@ -969,7 +969,7 @@ function ^(x::Float64, n::Integer) xnlo = ynlo = 0.0 if n < 0 rx = inv(x) - n==-2 && return x*x #keep compatability with literal_pow + n==-2 && return rx*rx #keep compatability with literal_pow isfinite(x) && (xnlo = -fma(x, rx, -1.) * rx) x = rx n = -n From 04cd0fb4372b2ba0987c2aea9234732f59013460 Mon Sep 17 00:00:00 2001 From: Oscar Smith Date: Mon, 8 Nov 2021 19:46:15 -0500 Subject: [PATCH 13/13] document reason for literal_pow change --- base/intfuncs.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/base/intfuncs.jl b/base/intfuncs.jl index 1bec3758b2334..9dc7e2ec67d0b 100644 --- a/base/intfuncs.jl +++ b/base/intfuncs.jl @@ -332,7 +332,7 @@ const HWNumber = Union{HWReal, Complex{<:HWReal}, Rational{<:HWReal}} @inline function literal_pow(f::typeof(^), x, ::Val{p}) where {p} if p < 0 if x isa BitInteger64 - f(Float64(x), p) + f(Float64(x), p) # inv would cause rounding, while Float64^Integer is able to compensate the inverse else f(inv(x), -p) end