Skip to content

Commit fe49d56

Browse files
authored
power uses Float64 exponents for integers (#53967)
Improve performance of `^(::Float64, n::Integer)` in the case of `abs(n) > 2^13`. While `pow_body` is unreliable for `abs(n) > 2^25` this implementation provides errors of a few ULPs, while runtime is capped to that of the `Float64` implementation. Fixes #53881 See also #53886.
1 parent c1b3661 commit fe49d56

File tree

4 files changed

+62
-13
lines changed

4 files changed

+62
-13
lines changed

base/math.jl

Lines changed: 41 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -1261,6 +1261,10 @@ function modf(x::T) where T<:IEEEFloat
12611261
return (rx, ix)
12621262
end
12631263

1264+
@inline function use_power_by_squaring(n::Integer)
1265+
-2^12 <= n <= 3 * 2^13
1266+
end
1267+
12641268
# @constprop aggressive to help the compiler see the switch between the integer and float
12651269
# variants for callers with constant `y`
12661270
@constprop :aggressive function ^(x::Float64, y::Float64)
@@ -1273,24 +1277,33 @@ end
12731277
y = sign(y)*0x1.8p62
12741278
end
12751279
yint = unsafe_trunc(Int64, y) # This is actually safe since julia freezes the result
1276-
y == yint && return @noinline x^yint
1277-
2*xu==0 && return abs(y)*Inf*(!(y>0)) # if x==0
1278-
x<0 && throw_exp_domainerror(x) # |y| is small enough that y isn't an integer
1279-
!isfinite(x) && return x*(y>0 || isnan(x)) # x is inf or NaN
1280+
yisint = y == yint
1281+
if yisint
1282+
yint == 0 && return 1.0
1283+
use_power_by_squaring(yint) && return @noinline pow_body(x, yint)
1284+
end
1285+
2*xu==0 && return abs(y)*Inf*(!(y>0)) # if x === +0.0 or -0.0 (Inf * false === 0.0)
1286+
s = 1
1287+
if x < 0
1288+
!yisint && throw_exp_domainerror(x) # y isn't an integer
1289+
s = ifelse(isodd(yint), -1, 1)
1290+
end
1291+
!isfinite(x) && return copysign(x,s)*(y>0 || isnan(x)) # x is inf or NaN
1292+
return copysign(pow_body(abs(x), y), s)
1293+
end
1294+
1295+
@assume_effects :foldable @noinline function pow_body(x::Float64, y::Float64)
1296+
xu = reinterpret(UInt64, x)
12801297
if xu < (UInt64(1)<<52) # x is subnormal
12811298
xu = reinterpret(UInt64, x * 0x1p52) # normalize x
12821299
xu &= ~sign_mask(Float64)
12831300
xu -= UInt64(52) << 52 # mess with the exponent
12841301
end
1285-
return pow_body(xu, y)
1286-
end
1287-
1288-
@inline function pow_body(xu::UInt64, y::Float64)
12891302
logxhi,logxlo = _log_ext(xu)
12901303
xyhi, xylo = two_mul(logxhi,y)
12911304
xylo = muladd(logxlo, y, xylo)
12921305
hi = xyhi+xylo
1293-
return Base.Math.exp_impl(hi, xylo-(hi-xyhi), Val(:ℯ))
1306+
return @inline Base.Math.exp_impl(hi, xylo-(hi-xyhi), Val(:ℯ))
12941307
end
12951308

12961309
@constprop :aggressive function ^(x::T, y::T) where T <: Union{Float16, Float32}
@@ -1314,12 +1327,29 @@ end
13141327
return T(exp2(log2(abs(widen(x))) * y))
13151328
end
13161329

1317-
# compensated power by squaring
13181330
@constprop :aggressive @inline function ^(x::Float64, n::Integer)
1331+
x^clamp(n, Int64)
1332+
end
1333+
@constprop :aggressive @inline function ^(x::Float64, n::Int64)
13191334
n == 0 && return one(x)
1320-
return pow_body(x, n)
1335+
if use_power_by_squaring(n)
1336+
return pow_body(x, n)
1337+
else
1338+
s = ifelse(x < 0 && isodd(n), -1.0, 1.0)
1339+
x = abs(x)
1340+
y = float(n)
1341+
if y == n
1342+
return copysign(pow_body(x, y), s)
1343+
else
1344+
n2 = n % 1024
1345+
y = float(n - n2)
1346+
return pow_body(x, y) * copysign(pow_body(x, n2), s)
1347+
end
1348+
end
13211349
end
13221350

1351+
# compensated power by squaring
1352+
# this method is only reliable for -2^20 < n < 2^20 (cf. #53881 #53886)
13231353
@assume_effects :terminates_locally @noinline function pow_body(x::Float64, n::Integer)
13241354
y = 1.0
13251355
xnlo = ynlo = 0.0

base/special/exp.jl

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -250,7 +250,7 @@ end
250250
twopk = (k + UInt64(53)) << 52
251251
return reinterpret(T, twopk + reinterpret(UInt64, small_part))*0x1p-53
252252
end
253-
#k == 1024 && return (small_part * 2.0) * 2.0^1023
253+
k == 1024 && return (small_part * 2.0) * 2.0^1023
254254
end
255255
twopk = Int64(k) << 52
256256
return reinterpret(T, twopk + reinterpret(Int64, small_part))

test/compiler/codegen.jl

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -866,7 +866,7 @@ if Sys.ARCH === :x86_64
866866
foo52079() = Core.Intrinsics.have_fma(Float64)
867867
if foo52079() == true
868868
let io = IOBuffer()
869-
code_native(io,^,(Float64,Float64), dump_module=false)
869+
code_native(io,Base.Math.exp_impl,(Float64,Float64,Val{:ℯ}), dump_module=false)
870870
str = String(take!(io))
871871
@test !occursin("fma_emulated", str)
872872
@test occursin("vfmadd", str)

test/math.jl

Lines changed: 19 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1470,6 +1470,25 @@ end
14701470
# two cases where we have observed > 1 ULP in the past
14711471
@test 0.0013653274095082324^-97.60372292227069 == 4.088393948750035e279
14721472
@test 8.758520413376658e-5^70.55863059215994 == 5.052076767078296e-287
1473+
1474+
# issue #53881
1475+
c53881 = 2.2844135865398217e222 # check correctness within 2 ULPs
1476+
@test prevfloat(1.0) ^ -Int64(2)^62 c53881 atol=2eps(c53881)
1477+
@test 2.0 ^ typemin(Int) == 0.0
1478+
@test (-1.0) ^ typemin(Int) == 1.0
1479+
Z = Int64(2)
1480+
E = prevfloat(1.0)
1481+
@test E ^ (-Z^54) 7.38905609893065
1482+
@test E ^ (-Z^62) 2.2844135865231613e222
1483+
@test E ^ (-Z^63) == Inf
1484+
@test abs(E ^ (Z^62-1) * E ^ (-Z^62+1) - 1) <= eps(1.0)
1485+
n, x = -1065564664, 0.9999997040311492
1486+
@test abs(x^n - Float64(big(x)^n)) / eps(x^n) == 0 # ULPs
1487+
@test E ^ (big(2)^100 + 1) == 0
1488+
@test E ^ 6705320061009595392 == nextfloat(0.0)
1489+
n = Int64(1024 / log2(E))
1490+
@test E^n == Inf
1491+
@test E^float(n) == Inf
14731492
end
14741493

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

0 commit comments

Comments
 (0)