don't inline float64^float64 (JuliaLang#42966)

Co-authored-by: Kristoffer Carlsson <kcarlsson89@gmail.com>

base/intfuncs.jl

@@ -319,6 +319,9 @@ const HWNumber = Union{HWReal, Complex{<:HWReal}, Rational{<:HWReal}}
 @inline literal_pow(::typeof(^), x::HWNumber, ::Val{1}) = x
 @inline literal_pow(::typeof(^), x::HWNumber, ::Val{2}) = x*x
 @inline literal_pow(::typeof(^), x::HWNumber, ::Val{3}) = x*x*x
+@inline literal_pow(::typeof(^), x::HWNumber, ::Val{-1}) = inv(x)
+@inline literal_pow(::typeof(^), x::HWNumber, ::Val{-2}) = (i=inv(x); i*i)
+@inline literal_pow(::typeof(^), x::HWNumber, ::Val{-3}) = (i=inv(x); i*i*i)
 
 # don't use the inv(x) transformation here since float^p is slightly more accurate
 @inline literal_pow(::typeof(^), x::AbstractFloat, ::Val{p}) where {p} = x^p
@@ -328,7 +331,11 @@ const HWNumber = Union{HWReal, Complex{<:HWReal}, Rational{<:HWReal}}
 # be computed in a type-stable way even for e.g. integers.
 @inline function literal_pow(f::typeof(^), x, ::Val{p}) where {p}
     if p < 0
-        literal_pow(^, inv(x), Val(-p))
+        if x isa BitInteger64
+            f(Float64(x), p) # inv would cause rounding, while Float64^Integer is able to compensate the inverse
+        else
+            f(inv(x), -p)
+        end
     else
         f(x, p)
     end

base/math.jl

@@ -939,7 +939,7 @@ function modf(x::T) where T<:IEEEFloat
     return (rx, ix)
 end
 
-@inline function ^(x::Float64, y::Float64)
+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
     x<0 && y > -4e18 && throw_exp_domainerror(x) # |y| is small enough that y isn't an integer
@@ -952,7 +952,7 @@ end
     hi = xyhi+xylo
     return Base.Math.exp_impl(hi, xylo-(hi-xyhi), Val(:ℯ))
 end
-@inline function ^(x::T, y::T) where T <: Union{Float16, Float32}
+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
     x < 0 && y > -4e18 && throw_exp_domainerror(x) # |y| is small enough that y isn't an integer
@@ -963,17 +963,18 @@ end
 end
 
 # compensated power by squaring
-@inline function ^(x::Float64, n::Integer)
+function ^(x::Float64, n::Integer)
     n == 0 && return one(x)
     y = 1.0
     xnlo = ynlo = 0.0
     if n < 0
         rx = inv(x)
+        n==-2 && return rx*rx #keep compatability with literal_pow
         isfinite(x) && (xnlo = -fma(x, rx, -1.) * rx)
         x = rx
         n = -n
     end
-    n==3 && return x*x*x #keep compatability with literal_pow
+    n == 3 && return x*x*x # keep compatibility with literal_pow
     while n > 1
         if n&1 > 0
             yn = x*y
@@ -988,9 +989,9 @@ end
     !isfinite(x) && return x*y
     return muladd(x, y, muladd(y, xnlo, x*ynlo))
 end
-@inline function ^(x::Float32, n::Integer)
+function ^(x::Float32, n::Integer)
     n < 0 && return inv(x)^(-n)
-    n==3 && return x*x*x #keep compatability with literal_pow
+    n == 3 && return x*x*x #keep compatibility with literal_pow
     Float32(Base.power_by_squaring(Float64(x),n))
 end
 @inline ^(x::Float16, y::Integer) = Float16(Float32(x) ^ y)

base/mathconstants.jl

@@ -76,7 +76,7 @@ julia> Base.MathConstants.eulergamma
 julia> dx = 10^-6;
 
 julia> sum(-exp(-x) * log(x) for x in dx:dx:100) * dx
-0.5772078382499134
+0.5772078382499133
 ```
 """
 γ, const eulergamma = γ