Skip to content

Commit 638a6a2

Browse files
fix Float32/Float16 power for giant integer exponents (#57829)
alternative to #57488 --------- Co-authored-by: Lilith Orion Hafner <[email protected]>
1 parent 3360a44 commit 638a6a2

File tree

2 files changed

+61
-11
lines changed

2 files changed

+61
-11
lines changed

base/math.jl

Lines changed: 33 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -1173,6 +1173,8 @@ end
11731173
return @inline Base.Math.exp_impl(hi, xylo-(hi-xyhi), Val(:ℯ))
11741174
end
11751175

1176+
# @constprop aggressive to help the compiler see the switch between the integer and float
1177+
# variants for callers with constant `y`
11761178
@constprop :aggressive function ^(x::T, y::T) where T <: Union{Float16, Float32}
11771179
x == 1 && return one(T)
11781180
# Exponents greater than this will always overflow or underflow.
@@ -1183,17 +1185,31 @@ end
11831185
y = sign(y)*max_exp
11841186
end
11851187
yint = unsafe_trunc(Int32, y) # This is actually safe since julia freezes the result
1186-
y == yint && return x^yint
1187-
x < 0 && throw_exp_domainerror(x)
1188-
!isfinite(x) && return x*(y>0 || isnan(x))
1189-
x==0 && return abs(y)*T(Inf)*(!(y>0))
1190-
return pow_body(x, y)
1188+
yisint = y == yint
1189+
if yisint
1190+
yint == 0 && return one(T)
1191+
use_power_by_squaring(yint) && return pow_body(x, yint)
1192+
end
1193+
s = 1
1194+
if x < 0
1195+
!yisint && throw_exp_domainerror(x) # y isn't an integer
1196+
s = ifelse(isodd(yint), -1, 1)
1197+
end
1198+
!isfinite(x) && return copysign(x,s)*(y>0 || isnan(x)) # x is inf or NaN
1199+
return copysign(pow_body(abs(x), y), s)
11911200
end
11921201

1193-
@inline function pow_body(x::T, y::T) where T <: Union{Float16, Float32}
1202+
@inline function pow_body(x::T, y) where T <: Union{Float16, Float32}
11941203
return T(exp2(log2(abs(widen(x))) * y))
11951204
end
11961205

1206+
@inline function pow_body(x::Union{Float16, Float32}, n::Int32)
1207+
n == -2 && return (i=inv(x); i*i)
1208+
n == 3 && return x*x*x #keep compatibility with literal_pow
1209+
n < 0 && return oftype(x, Base.power_by_squaring(inv(widen(x)), -n))
1210+
return oftype(x, Base.power_by_squaring(widen(x), n))
1211+
end
1212+
11971213
@constprop :aggressive @inline function ^(x::Float64, n::Integer)
11981214
n = clamp(n, Int64)
11991215
n == 0 && return one(x)
@@ -1242,11 +1258,17 @@ end
12421258
return ifelse(isfinite(x) & isfinite(err), muladd(x, y, err), x*y)
12431259
end
12441260

1245-
function ^(x::Union{Float16,Float32}, n::Integer)
1246-
n == -2 && return (i=inv(x); i*i)
1247-
n == 3 && return x*x*x #keep compatibility with literal_pow
1248-
n < 0 && return oftype(x, Base.power_by_squaring(inv(widen(x)),-n))
1249-
oftype(x, Base.power_by_squaring(widen(x),n))
1261+
# @constprop aggressive to help the compiler see the switch between the integer and float
1262+
# variants for callers with constant `y`
1263+
@constprop :aggressive @inline function ^(x::T, n::Integer) where T <: Union{Float16, Float32}
1264+
n = clamp(n, Int32)
1265+
# Exponents greater than this will always overflow or underflow.
1266+
# Note that NaN can pass through this, but that will end up fine.
1267+
n == 0 && return one(x)
1268+
use_power_by_squaring(n) && return pow_body(x, n)
1269+
s = ifelse(x < 0 && isodd(n), -one(T), one(T))
1270+
x = abs(x)
1271+
return pow_body(x, widen(T)(n))
12501272
end
12511273

12521274
## rem2pi-related calculations ##

test/math.jl

Lines changed: 28 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -1443,6 +1443,7 @@ end
14431443
Float32=>[.51, .51, .51, 2.0, 1.5],
14441444
Float64=>[.55, 0.8, 1.5, 2.0, 1.5])
14451445
for T in (Float16, Float32, Float64)
1446+
@inferred T T(1.1)^T(1.1) #test that we always return the right type
14461447
for x in (0.0, -0.0, 1.0, 10.0, 2.0, Inf, NaN, -Inf, -NaN)
14471448
for y in (0.0, -0.0, 1.0, -3.0,-10.0 , Inf, NaN, -Inf, -NaN)
14481449
got, expected = T(x)^T(y), T(big(x)^T(y))
@@ -1547,6 +1548,33 @@ end
15471548
@test EvenInteger(16) === @inferred EvenInteger(2)^Int(4) # avoid `literal_pow`
15481549
@test EvenInteger(16) === @inferred EvenInteger(2)^EvenInteger(4)
15491550
end
1551+
1552+
# issue #57464
1553+
@test Float32(1.1)^typemin(Int) == Float32(0.0)
1554+
@test Float16(1.1)^typemin(Int) == Float16(0.0)
1555+
@test Float32(1.1)^unsigned(0) === Float32(1.0)
1556+
@test Float32(1.1)^big(0) === Float32(1.0)
1557+
1558+
# By using a limited-precision integer (3 bits) we can trigger issue 57464
1559+
# for a case where the answer isn't zero.
1560+
struct Int3 <: Integer
1561+
x::Int8
1562+
function Int3(x::Integer)
1563+
if x < -4 || x > 3
1564+
Core.throw_inexacterror(:Int3, Int3, x)
1565+
end
1566+
return new(x)
1567+
end
1568+
end
1569+
Base.typemin(::Type{Int3}) = Int3(-4)
1570+
Base.promote_rule(::Type{Int3}, ::Type{T}) where {T<:Integer} = T
1571+
Base.convert(::Type{T}, x::Int3) where {T<:Integer} = convert(T, x.x)
1572+
Base.:-(x::Int3) = x.x == -4 ? x : Int3(-x.x)
1573+
Base.trailing_zeros(x::Int3) = trailing_zeros(x.x)
1574+
Base.:>>(x::Int3, n::UInt64) = Int3(x.x>>n)
1575+
1576+
@test 1.001f0^-3 == 1.001f0^Int3(-3)
1577+
@test 1.001f0^-4 == 1.001f0^typemin(Int3)
15501578
end
15511579

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

0 commit comments

Comments
 (0)