fix precision issue in Float64^Float64.

base/intfuncs.jl

@@ -321,7 +321,6 @@ const HWNumber = Union{HWReal, Complex{<:HWReal}, Rational{<:HWReal}}
 @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

base/math.jl

@@ -1010,8 +1010,8 @@ end
     !isfinite(x) && return x*(y>0 || isnan(x))
     x==0 && return abs(y)*Inf*(!(y>0))
     logxhi,logxlo = Base.Math._log_ext(x)
-    xyhi = logxhi*y
-    xylo = logxlo*y
+    xyhi, xylo = two_mul(logxhi,y)
+    xylo = muladd(logxlo, y, xylo)
     hi = xyhi+xylo
     return Base.Math.exp_impl(hi, xylo-(hi-xyhi), Val(:ℯ))
 end
@@ -1041,14 +1041,14 @@ end
 @assume_effects :terminates_locally @noinline function pow_body(x::Float64, n::Integer)
     y = 1.0
     xnlo = ynlo = 0.0
+    n == 3 && return x*x*x # keep compatibility with literal_pow
     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 compatibility with literal_pow
     while n > 1
         if n&1 > 0
             err = muladd(y, xnlo, x*ynlo)
@@ -1065,8 +1065,9 @@ end
 end
 
 function ^(x::Float32, n::Integer)
-    n < 0 && return inv(x)^(-n)
+    n == -2 && return (i=inv(x); i*i)
     n == 3 && return x*x*x #keep compatibility with literal_pow
+    n < 0 && return Float32(Base.power_by_squaring(inv(Float64(x)),-n))
     Float32(Base.power_by_squaring(Float64(x),n))
 end
 @inline ^(x::Float16, y::Integer) = Float16(Float32(x) ^ y)

doc/src/manual/complex-and-rational-numbers.md

@@ -36,7 +36,7 @@ julia> (-1 + 2im)^2
 -3 - 4im
 
 julia> (-1 + 2im)^2.5
-2.7296244647840084 - 6.960664459571898im
+2.729624464784009 - 6.9606644595719im
 
 julia> (-1 + 2im)^(1 + 1im)
 -0.27910381075826657 + 0.08708053414102428im

test/math.jl

@@ -155,12 +155,6 @@ 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
             @test (-x)^(yi+1) == -(x^(yi+1))
             @test acos(x) ≈ acos(big(x))
@@ -1323,6 +1317,30 @@ end
     end
 end
 
+@testset "pow" begin
+    for T in (Float16, Float32, Float64)
+        for x in (0.0, -0.0, 1.0, 10.0, 2.0, Inf, NaN, -Inf, -NaN)
+            for y in (0.0, -0.0, 1.0, -3.0,-10.0 , Inf, NaN, -Inf, -NaN)
+                got, expected = T(x)^T(y), T(big(x))^T(y)
+                @test isnan_type(T, got) && isnan_type(T, expected) || (got === expected)
+            end
+        end
+        for _ in 1:2^16
+            x=rand(T)*100; y=rand(T)*200-100
+            got, expected = x^y, widen(x)^y
+            if isfinite(eps(T(expected)))
+                @test abs(expected-got) <= 1.3*eps(T(expected)) || (x,y)
+            end
+        end
+        # 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)
+    end
+    # test for large negative exponent where error compensation matters
+    @test 0.9999999955206014^-1.0e8 == 1.565084574870928
+end
+
 # Test that sqrt behaves correctly and doesn't exhibit fp80 double rounding.
 # This happened on old glibc versions.
 # Test case from https://sourceware.org/bugzilla/show_bug.cgi?id=14032.