fix Float32/Float16 power for giant integer exponents

base/math.jl

@@ -1134,22 +1134,21 @@ 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)
-    xu = reinterpret(UInt64, x)
-    xu == reinterpret(UInt64, 1.0) && return 1.0
+@constprop :aggressive function ^(x::T, y::T) where T <: Union{Float16, Float32}
+    x == 1 && return one(T)
     # Exponents greater than this will always overflow or underflow.
     # Note that NaN can pass through this, but that will end up fine.
-    if !(abs(y)<0x1.8p62)
+    max_exp = T == Float16 ? T(3<<14) : T(0x1.Ap30)
+    if !(abs(y)<max_exp)
         isnan(y) && return y
-        y = sign(y)*0x1.8p62
+        y = sign(y)*max_exp
     end
-    yint = unsafe_trunc(Int64, y) # This is actually safe since julia freezes the result
+    yint = unsafe_trunc(Int32, y) # This is actually safe since julia freezes the result
     yisint = y == yint
     if yisint
         yint == 0 && return 1.0
-        use_power_by_squaring(yint) && return @noinline pow_body(x, yint)
+        use_power_by_squaring(yint) && return 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
@@ -1159,39 +1158,55 @@ end
     return copysign(pow_body(abs(x), y), s)
 end
 
-@assume_effects :foldable @noinline function pow_body(x::Float64, y::Float64)
-    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
-    logxhi,logxlo = _log_ext(xu)
-    xyhi, xylo = two_mul(logxhi,y)
-    xylo = muladd(logxlo, y, xylo)
-    hi = xyhi+xylo
-    return @inline Base.Math.exp_impl(hi, xylo-(hi-xyhi), Val(:ℯ))
+# @constprop aggressive to help the compiler see the switch between the integer and float
+# variants for callers with constant `y`
+@constprop :aggressive @inline function ^(x::T, n::Integer) where T <: Union{Float16, Float32}
+    n = clamp(n, Int32)
+    # Exponents greater than this will always overflow or underflow.
+    # Note that NaN can pass through this, but that will end up fine.
+    n == 0 && return one(x)
+    use_power_by_squaring(n) && return pow_body(x, n)
+    s = ifelse(x < 0 && isodd(n), -1f0, 1f0)
+    x = abs(x)
+    return pow_body(x, widen(T)(n))
 end
 
-@constprop :aggressive function ^(x::T, y::T) where T <: Union{Float16, Float32}
-    x == 1 && return one(T)
+@inline function pow_body(x::T, y) where T <: Union{Float16, Float32}
+    return T(exp2(log2(abs(widen(x))) * y))
+end
+
+@inline function pow_body(x::Union{Float16, Float32}, n::Int32)
+    n == -2 && return (i=inv(x); i*i)
+    n == 3 && return x*x*x #keep compatibility with literal_pow
+    n < 0 && return oftype(x, Base.power_by_squaring(inv(widen(x)), -n))
+    return oftype(x, Base.power_by_squaring(widen(x), n))
+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)
+    xu = reinterpret(UInt64, x)
+    xu == reinterpret(UInt64, 1.0) && return 1.0
     # Exponents greater than this will always overflow or underflow.
     # Note that NaN can pass through this, but that will end up fine.
-    max_exp = T == Float16 ? T(3<<14) : T(0x1.Ap30)
-    if !(abs(y)<max_exp)
+    if !(abs(y)<0x1.8p62)
         isnan(y) && return y
-        y = sign(y)*max_exp
+        y = sign(y)*0x1.8p62
     end
-    yint = unsafe_trunc(Int32, y) # This is actually safe since julia freezes the result
-    y == yint && return x^yint
-    x < 0 && throw_exp_domainerror(x)
-    !isfinite(x) && return x*(y>0 || isnan(x))
-    x==0 && return abs(y)*T(Inf)*(!(y>0))
-    return pow_body(x, y)
-end
-
-@inline function pow_body(x::T, y::T) where T <: Union{Float16, Float32}
-    return T(exp2(log2(abs(widen(x))) * y))
+    yint = unsafe_trunc(Int64, y) # This is actually safe since julia freezes the result
+    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)
 end
 
 @constprop :aggressive @inline function ^(x::Float64, n::Integer)
@@ -1213,6 +1228,20 @@ end
     end
 end
 
+@assume_effects :foldable @noinline function pow_body(x::Float64, y::Float64)
+    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
+    logxhi,logxlo = _log_ext(xu)
+    xyhi, xylo = two_mul(logxhi,y)
+    xylo = muladd(logxlo, y, xylo)
+    hi = xyhi+xylo
+    return @inline Base.Math.exp_impl(hi, xylo-(hi-xyhi), Val(:ℯ))
+end
+
 # compensated power by squaring
 # 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)
@@ -1242,13 +1271,6 @@ end
     return ifelse(isfinite(x) & isfinite(err), muladd(x, y, err), x*y)
 end
 
-function ^(x::Union{Float16,Float32}, n::Integer)
-    n == -2 && return (i=inv(x); i*i)
-    n == 3 && return x*x*x #keep compatibility with literal_pow
-    n < 0 && return oftype(x, Base.power_by_squaring(inv(widen(x)),-n))
-    oftype(x, Base.power_by_squaring(widen(x),n))
-end
-
 ## rem2pi-related calculations ##
 
 function add22condh(xh::Float64, xl::Float64, yh::Float64, yl::Float64)

test/math.jl

@@ -1547,6 +1547,33 @@ end
         @test EvenInteger(16) === @inferred EvenInteger(2)^Int(4)  # avoid `literal_pow`
         @test EvenInteger(16) === @inferred EvenInteger(2)^EvenInteger(4)
     end
+
+    # issue #57464
+    @test Float32(1.1)^typemin(Int) == Float32(0.0)
+    @test Float16(1.1)^typemin(Int) == Float16(0.0)
+    @test Float32(1.1)^unsigned(0) === Float32(1.0)
+    @test Float32(1.1)^big(0) === Float32(1.0)
+
+    # By using a limited-precision integer (3 bits) we can trigger issue 57464
+    # for a case where the answer isn't zero.
+    struct Int3 <: Integer
+        x::Int8
+        function Int3(x::Integer)
+            if x < -4 || x > 3
+                Core.throw_inexacterror(:Int3, Int3, x)
+            end
+            return new(x)
+        end
+    end
+    Base.typemin(::Type{Int3}) = Int3(-4)
+    Base.promote_rule(::Type{Int3}, ::Type{Int}) = Int
+    Base.convert(::Type{Int64}, x::Int3) = convert(Int64, x.x)
+    Base.:-(x::Int3) = x.x == -4 ? x : Int3(-x.x)
+    Base.trailing_zeros(x::Int3) = trailing_zeros(x.x)
+    Base.:>>(x::Int3, n::UInt64) = Int3(x.x>>n)
+
+    @test 1.001f0^-3 == 1.001f0^Int3(-3)
+    @test 1.001f0^-4 == 1.001f0^typemin(Int3)
 end
 
 @testset "special function `::Real` fallback shouldn't recur without bound, issue #57789" begin