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.

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

power uses Float64 exponents for integers #53967

Merged
merged 25 commits into from
Apr 22, 2024
Merged
Show file tree
Hide file tree
Changes from 8 commits
Commits
Show all changes
25 commits
Select commit Hold shift + click to select a range
81a0813
power uses Float64 exponents for integers
KlausC Apr 5, 2024
5c60c57
use_pbs - better runtime prediction
KlausC Apr 6, 2024
c861bb0
fixed lead_zeros bug introduced previously
KlausC Apr 6, 2024
4ab7026
fixed corner cases and added tests
KlausC Apr 10, 2024
bb7633c
use_pbs threshold
KlausC Apr 10, 2024
fed55b7
fixed effects in `^`
KlausC Apr 11, 2024
fd7fc04
fix tests for pow
KlausC Apr 12, 2024
bfbd7d3
Merge branch 'master' into krc/pow-float-path
KlausC Apr 12, 2024
9f8e198
lest @assume_effects
KlausC Apr 12, 2024
703c047
Merge branch 'master' into krc/pow-float-path
KlausC Apr 13, 2024
ffe1c25
simplified `use_power_by_squaring`
KlausC Apr 13, 2024
df1091d
Merge branch 'krc/pow-float-path' of https://github.com/KlausC/julia …
KlausC Apr 13, 2024
e9b33a2
adjusted nmin, nmax
KlausC Apr 13, 2024
06d0a7c
use `clamp`
KlausC Apr 13, 2024
7212b2d
Merge branch 'master' into krc/pow-float-path
KlausC Apr 13, 2024
09658fc
Merge branch 'master' into krc/pow-float-path
KlausC Apr 14, 2024
0964a16
Merge branch 'master' into krc/pow-float-path
KlausC Apr 15, 2024
e43797d
Merge branch 'master' into krc/pow-float-path
KlausC Apr 16, 2024
d40a3c8
Merge branch 'master' into krc/pow-float-path
KlausC Apr 17, 2024
9f6d39d
try to inline `Base.Math.exp_impl`
KlausC Apr 18, 2024
f6c24e1
Merge branch 'master' into krc/pow-float-path
KlausC Apr 18, 2024
2d6fdf5
noinline pow_body, other compiler test
KlausC Apr 18, 2024
249c53f
Merge branch 'master' into krc/pow-float-path
KlausC Apr 19, 2024
2260a80
Merge branch 'master' into krc/pow-float-path
KlausC Apr 20, 2024
4c7f788
removed redundant code, mitigate test case
KlausC Apr 21, 2024
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
63 changes: 52 additions & 11 deletions base/math.jl
Original file line number Diff line number Diff line change
Expand Up @@ -37,7 +37,7 @@ end
throw(DomainError(x,
LazyString(f," was called with a real argument < -1 but will only return a complex result if called with a complex argument. Try ", f,"(Complex(x)).")))
end
@noinline function throw_exp_domainerror(x)
@assume_effects :terminates_globally @noinline function throw_exp_domainerror(x)
KlausC marked this conversation as resolved.
Show resolved Hide resolved
throw(DomainError(x, LazyString(
KlausC marked this conversation as resolved.
Show resolved Hide resolved
"Exponentiation yielding a complex result requires a ",
"complex argument.\nReplace x^y with (x+0im)^y, ",
Expand Down Expand Up @@ -1255,6 +1255,15 @@ function modf(x::T) where T<:IEEEFloat
return (rx, ix)
end

@inline function use_power_by_squaring(n::Union{Int32,Int64,Int128,UInt32,UInt64,UInt128})
KlausC marked this conversation as resolved.
Show resolved Hide resolved
# top_set_bit is not available during bootstrap
x = abs(n)
2*(8sizeof(x) - leading_zeros(x)) + count_ones(x) + 6*(n<0) < 40
end
@inline function use_power_by_squaring(n::Integer)
-2^13 <= n < 2^13
end

# @constprop aggressive to help the compiler see the switch between the integer and float
# variants for callers with constant `y`
@constprop :aggressive function ^(x::Float64, y::Float64)
Expand All @@ -1267,19 +1276,28 @@ end
y = sign(y)*0x1.8p62
end
yint = unsafe_trunc(Int64, y) # This is actually safe since julia freezes the result
y == yint && return @noinline x^yint
2*xu==0 && return abs(y)*Inf*(!(y>0)) # if x==0
x<0 && throw_exp_domainerror(x) # |y| is small enough that y isn't an integer
!isfinite(x) && return x*(y>0 || isnan(x)) # x is inf or NaN
yisint = y == yint
if yisint
yint == 0 && return 1.0
use_power_by_squaring(yint) && return @noinline pow_body(x, yint)
end
2*xu==0 && return abs(y)*Inf*(!(y>0)) # if x === +0.0 or -0.0 (Inf * false === 0.0)
s = 1
if x < 0
!yisint && throw_exp_domainerror(x) # y isn't an integer
s = ifelse(isodd(yint), -1, 1)
end
!isfinite(x) && return copysign(x,s)*(y>0 || isnan(x)) # x is inf or NaN
return copysign(pow_body(abs(x), y), s)
oscardssmith marked this conversation as resolved.
Show resolved Hide resolved
end

@assume_effects :total @noinline function pow_body(x::Float64, y::Float64)
KlausC marked this conversation as resolved.
Show resolved Hide resolved
xu = reinterpret(UInt64, x)
if xu < (UInt64(1)<<52) # x is subnormal
xu = reinterpret(UInt64, x * 0x1p52) # normalize x
xu &= ~sign_mask(Float64)
xu -= UInt64(52) << 52 # mess with the exponent
end
return pow_body(xu, y)
end

@inline function pow_body(xu::UInt64, y::Float64)
logxhi,logxlo = _log_ext(xu)
xyhi, xylo = two_mul(logxhi,y)
xylo = muladd(logxlo, y, xylo)
Expand Down Expand Up @@ -1308,14 +1326,37 @@ end
return T(exp2(log2(abs(widen(x))) * y))
end

# compensated power by squaring
@constprop :aggressive @inline function ^(x::Float64, n::Integer)
T = Int64
m = n >= typemax(T) ? typemax(T) : n <= typemin(T) ? typemin(T) : n
x^(m%Int64)
end
KlausC marked this conversation as resolved.
Show resolved Hide resolved
# compensated power by squaring
@constprop :aggressive @inline function ^(x::Float64, n::Int64)
n == 0 && return one(x)
return pow_body(x, n)
if use_power_by_squaring(n)
return pow_body(x, n)
else
s = ifelse(x < 0 && isodd(n), -1.0, 1.0)
KlausC marked this conversation as resolved.
Show resolved Hide resolved
x = abs(x)
y = Float64(n)
if y == n
return copysign(pow_body(x, y), s)
else
if n % Int64 != n
n = ifelse(n < 0, typemin(Int64), typemax(Int64))
end
n2 = n % 1024
y = float(n - n2)
return pow_body(x, y) * copysign(pow_body(x, n2), s)
end
end
end

# this method is only reliable for -2^20 < n < 2^20 (cf. #53881 #53886)
@assume_effects :terminates_locally @noinline function pow_body(x::Float64, n::Integer)
y = 1.0
n == 0 && return y
KlausC marked this conversation as resolved.
Show resolved Hide resolved
xnlo = ynlo = 0.0
n == 3 && return x*x*x # keep compatibility with literal_pow
if n < 0
Expand Down
2 changes: 1 addition & 1 deletion base/special/exp.jl
Original file line number Diff line number Diff line change
Expand Up @@ -250,7 +250,7 @@ end
twopk = (k + UInt64(53)) << 52
return reinterpret(T, twopk + reinterpret(UInt64, small_part))*0x1p-53
end
#k == 1024 && return (small_part * 2.0) * 2.0^1023
k == 1024 && return (small_part * 2.0) * 2.0^1023
oscardssmith marked this conversation as resolved.
Show resolved Hide resolved
end
twopk = Int64(k) << 52
return reinterpret(T, twopk + reinterpret(Int64, small_part))
Expand Down
18 changes: 18 additions & 0 deletions test/math.jl
Original file line number Diff line number Diff line change
Expand Up @@ -1459,6 +1459,24 @@ 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

# issue #53881
@test prevfloat(1.0) ^ -Int64(2)^62 == 2.2844135865398217e222
KlausC marked this conversation as resolved.
Show resolved Hide resolved
@test 2.0 ^ typemin(Int) == 0.0
@test (-1.0) ^ typemin(Int) == 1.0
Z = Int64(2)
E = prevfloat(1.0)
@test E ^ (-Z^54) ≈ 7.38905609893065
@test E ^ (-Z^62) ≈ 2.2844135865231613e222
@test E ^ (-Z^63) == Inf
@test abs(E ^ (Z^62-1) * E ^ (-Z^62+1) - 1) <= eps(1.0)
n, x = -1065564664, 0.9999997040311492
@test abs(x^n - Float64(big(x)^n)) / eps(x^n) == 0 # ULPs
@test E ^ (big(2)^100 + 1) == 0
@test E ^ 6705320061009595392 == nextfloat(0.0)
n = Int64(1024 / log2(E))
@test E^n == Inf
@test E^float(n) == Inf
end

# Test that sqrt behaves correctly and doesn't exhibit fp80 double rounding.
Expand Down