Skip to content

Commit cbc56b0

Browse files
committed
Revert "power uses Float64 exponents for integers (#53967)"
This reverts commit fe073c4.
1 parent 95c2125 commit cbc56b0

File tree

3 files changed

+13
-62
lines changed

3 files changed

+13
-62
lines changed

base/math.jl

Lines changed: 12 additions & 42 deletions
Original file line numberDiff line numberDiff line change
@@ -1189,10 +1189,6 @@ function modf(x::T) where T<:IEEEFloat
11891189
return (rx, ix)
11901190
end
11911191

1192-
@inline function use_power_by_squaring(n::Integer)
1193-
-2^12 <= n <= 3 * 2^13
1194-
end
1195-
11961192
# @constprop aggressive to help the compiler see the switch between the integer and float
11971193
# variants for callers with constant `y`
11981194
@constprop :aggressive function ^(x::Float64, y::Float64)
@@ -1205,33 +1201,24 @@ end
12051201
y = sign(y)*0x1.8p62
12061202
end
12071203
yint = unsafe_trunc(Int64, y) # This is actually safe since julia freezes the result
1208-
yisint = y == yint
1209-
if yisint
1210-
yint == 0 && return 1.0
1211-
use_power_by_squaring(yint) && return @noinline pow_body(x, yint)
1212-
end
1213-
2*xu==0 && return abs(y)*Inf*(!(y>0)) # if x === +0.0 or -0.0 (Inf * false === 0.0)
1214-
s = 1
1215-
if x < 0
1216-
!yisint && throw_exp_domainerror(x) # y isn't an integer
1217-
s = ifelse(isodd(yint), -1, 1)
1218-
end
1219-
!isfinite(x) && return copysign(x,s)*(y>0 || isnan(x)) # x is inf or NaN
1220-
return copysign(pow_body(abs(x), y), s)
1221-
end
1222-
1223-
@assume_effects :foldable @noinline function pow_body(x::Float64, y::Float64)
1224-
xu = reinterpret(UInt64, x)
1204+
y == yint && return @noinline x^yint
1205+
2*xu==0 && return abs(y)*Inf*(!(y>0)) # if x==0
1206+
x<0 && throw_exp_domainerror(x) # |y| is small enough that y isn't an integer
1207+
!isfinite(x) && return x*(y>0 || isnan(x)) # x is inf or NaN
12251208
if xu < (UInt64(1)<<52) # x is subnormal
12261209
xu = reinterpret(UInt64, x * 0x1p52) # normalize x
12271210
xu &= ~sign_mask(Float64)
12281211
xu -= UInt64(52) << 52 # mess with the exponent
12291212
end
1230-
logxhi,logxlo = _log_ext(xu)
1213+
return pow_body(xu, y)
1214+
end
1215+
1216+
@inline function pow_body(xu::UInt64, y::Float64)
1217+
logxhi,logxlo = Base.Math._log_ext(xu)
12311218
xyhi, xylo = two_mul(logxhi,y)
12321219
xylo = muladd(logxlo, y, xylo)
12331220
hi = xyhi+xylo
1234-
return @inline Base.Math.exp_impl(hi, xylo-(hi-xyhi), Val(:ℯ))
1221+
return Base.Math.exp_impl(hi, xylo-(hi-xyhi), Val(:ℯ))
12351222
end
12361223

12371224
@constprop :aggressive function ^(x::T, y::T) where T <: Union{Float16, Float32}
@@ -1255,29 +1242,12 @@ end
12551242
return T(exp2(log2(abs(widen(x))) * y))
12561243
end
12571244

1245+
# compensated power by squaring
12581246
@constprop :aggressive @inline function ^(x::Float64, n::Integer)
1259-
x^clamp(n, Int64)
1260-
end
1261-
@constprop :aggressive @inline function ^(x::Float64, n::Int64)
12621247
n == 0 && return one(x)
1263-
if use_power_by_squaring(n)
1264-
return pow_body(x, n)
1265-
else
1266-
s = ifelse(x < 0 && isodd(n), -1.0, 1.0)
1267-
x = abs(x)
1268-
y = float(n)
1269-
if y == n
1270-
return copysign(pow_body(x, y), s)
1271-
else
1272-
n2 = n % 1024
1273-
y = float(n - n2)
1274-
return pow_body(x, y) * copysign(pow_body(x, n2), s)
1275-
end
1276-
end
1248+
return pow_body(x, n)
12771249
end
12781250

1279-
# compensated power by squaring
1280-
# this method is only reliable for -2^20 < n < 2^20 (cf. #53881 #53886)
12811251
@assume_effects :terminates_locally @noinline function pow_body(x::Float64, n::Integer)
12821252
y = 1.0
12831253
xnlo = ynlo = 0.0

base/special/exp.jl

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -252,7 +252,7 @@ end
252252
twopk = (k + UInt64(53)) << 52
253253
return reinterpret(T, twopk + reinterpret(UInt64, small_part))*0x1p-53
254254
end
255-
k == 1024 && return (small_part * 2.0) * 0x1p1023
255+
# k == 1024 && return (small_part * 2.0) * 0x1p1023
256256
end
257257
twopk = Int64(k) << 52
258258
return reinterpret(T, twopk + reinterpret(Int64, small_part))

test/math.jl

Lines changed: 0 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -1453,25 +1453,6 @@ end
14531453
# two cases where we have observed > 1 ULP in the past
14541454
@test 0.0013653274095082324^-97.60372292227069 == 4.088393948750035e279
14551455
@test 8.758520413376658e-5^70.55863059215994 == 5.052076767078296e-287
1456-
1457-
# issue #53881
1458-
c53881 = 2.2844135865398217e222 # check correctness within 2 ULPs
1459-
@test prevfloat(1.0) ^ -Int64(2)^62 c53881 atol=2eps(c53881)
1460-
@test 2.0 ^ typemin(Int) == 0.0
1461-
@test (-1.0) ^ typemin(Int) == 1.0
1462-
Z = Int64(2)
1463-
E = prevfloat(1.0)
1464-
@test E ^ (-Z^54) 7.38905609893065
1465-
@test E ^ (-Z^62) 2.2844135865231613e222
1466-
@test E ^ (-Z^63) == Inf
1467-
@test abs(E ^ (Z^62-1) * E ^ (-Z^62+1) - 1) <= eps(1.0)
1468-
n, x = -1065564664, 0.9999997040311492
1469-
@test abs(x^n - Float64(big(x)^n)) / eps(x^n) == 0 # ULPs
1470-
@test E ^ (big(2)^100 + 1) == 0
1471-
@test E ^ 6705320061009595392 == nextfloat(0.0)
1472-
n = Int64(1024 / log2(E))
1473-
@test E^n == Inf
1474-
@test E^float(n) == Inf
14751456
end
14761457

14771458
@testset "special function `::Real` fallback shouldn't recur without bound, issue #57789" begin

0 commit comments

Comments
 (0)