Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

don't inline float64^float64 #42966

Merged
merged 13 commits into from
Nov 10, 2021
Merged

don't inline float64^float64 #42966

Show file tree
Hide file tree
Changes from 12 commits
Commits
Show all changes
13 commits
Select commit Hold shift + click to select a range
07ded3c
don't inline float64^float64
oscardssmith Nov 5, 2021
10c4a6b
fix literal_pow to be more accurate for large negative bases
oscardssmith Nov 8, 2021
bcf5e27
fixup! fix literal_pow to be more accurate for large negative bases
oscardssmith Nov 8, 2021
9a3c9db
fix matrix^int
oscardssmith Nov 8, 2021
0e7e160
restrict special case to small ints
oscardssmith Nov 8, 2021
e34a79d
fix broken
oscardssmith Nov 8, 2021
e3704bb
fix part 2
oscardssmith Nov 8, 2021
1ef574f
Update base/intfuncs.jl
oscardssmith Nov 8, 2021
be196d7
why is this so annoying
oscardssmith Nov 8, 2021
ac3677f
remove the last inlines.
KristofferC Nov 8, 2021
ce5b5f3
fix another bug
oscardssmith Nov 8, 2021
0d5a68c
inv(x)!=x
oscardssmith Nov 8, 2021
04cd0fb
document reason for literal_pow change
oscardssmith Nov 9, 2021
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
9 changes: 8 additions & 1 deletion base/intfuncs.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -328,7 +331,11 @@ 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
literal_pow(^, inv(x), Val(-p))
if x isa BitInteger64
f(Float64(x), p)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Maybe add some comment on this? Is this to improve performance, accuracy, or both?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Good idea. The reasoning is that computing inv(x)^(-p) rounds twice. Once to compute inv(x) and once to compute the power. Since the Float64^Float64 code works to half an ULP (even for negative numbers), we want to let that code handle the negative powers since it does so better.

else
f(inv(x), -p)
end
else
f(x, p)
end
Expand Down
13 changes: 7 additions & 6 deletions base/math.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand All @@ -952,7 +952,7 @@ end
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
Expand All @@ -963,17 +963,18 @@ 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
if n < 0
rx = inv(x)
n==-2 && return rx*rx #keep compatability with literal_pow
isfinite(x) && (xnlo = -fma(x, rx, -1.) * rx)
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
Expand All @@ -988,9 +989,9 @@ end
!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)
Expand Down
2 changes: 1 addition & 1 deletion base/mathconstants.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -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 = γ
Expand Down