fix negative numbers to powers >2^64 (#44456)

*fix negative numbers to powers >2^64

Author: oscardssmith

base/math.jl

@@ -999,6 +999,8 @@ end
 @constprop :aggressive function ^(x::Float64, y::Float64)
     yint = unsafe_trunc(Int, y) # Note, this is actually safe since julia freezes the result
     y == yint && return x^yint
+    #numbers greater than 2*inv(eps(T)) must be even, and the pow will overflow
+    y >= 2*inv(eps()) && return x^(typemax(Int64)-1)
     x<0 && y > -4e18 && throw_exp_domainerror(x) # |y| is small enough that y isn't an integer
     x == 1 && return 1.0
     return pow_body(x, y)
@@ -1017,6 +1019,8 @@ end
 @constprop :aggressive function ^(x::T, y::T) where T <: Union{Float16, Float32}
     yint = unsafe_trunc(Int64, y) # Note, this is actually safe since julia freezes the result
     y == yint && return x^yint
+    #numbers greater than 2*inv(eps(T)) must be even, and the pow will overflow
+    y >= 2*inv(eps(T)) && return x^(typemax(Int64)-1)
     x < 0 && y > -4e18 && throw_exp_domainerror(x) # |y| is small enough that y isn't an integer
     return pow_body(x, y)
 end

test/math.jl

@@ -155,6 +155,10 @@ end
             @test x^y === T(big(x)^big(y))
             @test x^1 === x
             @test x^yi === T(big(x)^yi)
+            # test (-x)^y for y larger than typemax(Int)
+            @test T(-1)^floatmax(T) === T(1)
+            @test prevfloat(T(-1))^floatmax(T) === T(Inf)
+            @test nextfloat(T(-1))^floatmax(T) === T(0.0)
             # test for large negative exponent where error compensation matters
             @test 0.9999999955206014^-1.0e8 == 1.565084574870928
             @test (-x)^yi == x^yi