diff --git a/base/math.jl b/base/math.jl
index c3c15505c5f8f..197f50ebd5946 100644
--- a/base/math.jl
+++ b/base/math.jl
@@ -1314,32 +1314,69 @@ end
     return pow_body(x, n)
 end
 
+# high accuracy implementation of power by squaring
+# calculations are effectively in double precision using fma
 @assume_effects :terminates_locally @noinline function pow_body(x::Float64, n::Integer)
     y = 1.0
-    xnlo = ynlo = 0.0
+    isnan(x) && return x
     n == 3 && return x*x*x # keep compatibility with literal_pow
+    normal(x, y) = begin z = x; x = z + y; y -= x - z; return x, y; end
+    xnlo = ynlo = 0.0
+    negate = signbit(x) && n & 1 > 0
+    xa = abs(x)
+    toinf = xa >= 1
+    invert = false
     if n < 0
-        rx = inv(x)
-        n==-2 && return rx*rx #keep compatibility with literal_pow
-        isfinite(x) && (xnlo = -fma(x, rx, -1.) * rx)
-        x = rx
+        if n == -2
+            rx = inv(x)
+            return rx * rx #keep compatibility with literal_pow
+        end
+        if toinf # || xa < 1.0 - sqrt(eps(1.0))
+            xx, x = x, inv(x)
+            xnlo = fma(-xx, x, 1.0) * x
+        else
+            invert = true
+        end
         n = -n
+        toinf = !toinf
     end
-    while n > 1
-        if n&1 > 0
-            err = muladd(y, xnlo, x*ynlo)
-            y, ynlo = two_mul(x,y)
-            ynlo += err
+    MAXI = 16
+    maxi = invert ? 1 : MAXI
+    while n != 0
+        x, xnlo = normal(x, xnlo)
+        x, xnlo, y, ynlo, n = pow_loop(x, xnlo, y, ynlo, n, maxi)
+
+        if invert
+            xx, x = x, inv(x)
+            xnlo = (fma(-xx, x, 1.0) - xnlo * x) * x
+            yy, y = y, inv(y)
+            ynlo = (fma(-yy, y, 1.0) - ynlo * y) * y
+            invert = false
+            maxi = MAXI
+        end
+        y, ynlo = normal(y, ynlo)
+    end
+    return isfinite(y) && !iszero(y) ? y : flipsign(ifelse(toinf, Inf, 0.0), -negate)
+end
+
+@assume_effects :terminates_locally @inline function pow_loop(x, xnlo, y, ynlo, n, maxi)
+    for i = 1:maxi
+        if n & 1 > 0
+            err = y * xnlo + x * ynlo
+            xx = y
+            y = x * xx
+            ynlo = fma(x, xx, -y) + err
         end
-        err = x*2*xnlo
-        x, xnlo = two_mul(x, x)
-        xnlo += err
         n >>>= 1
+        n == 0 && break
+        xx = x * x
+        xnlo = fma(x, x, -xx) + x * 2 * xnlo
+        x = xx
     end
-    err = muladd(y, xnlo, x*ynlo)
-    return ifelse(isfinite(x) & isfinite(err), muladd(x, y, err), x*y)
+    x, xnlo, y, ynlo, n
 end
 
+
 function ^(x::Float32, n::Integer)
     n == -2 && return (i=inv(x); i*i)
     n == 3 && return x*x*x #keep compatibility with literal_pow
diff --git a/test/math.jl b/test/math.jl
index f551bb3a5d4b7..179149e321996 100644
--- a/test/math.jl
+++ b/test/math.jl
@@ -1454,11 +1454,22 @@ end
         @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
+    @test prevfloat(1.0) ^ -Int64(2)^62 == 2.2844135865398217e222
     @test 3e18^20 == Inf
     # two cases where we have observed > 1 ULP in the past
     @test 0.0013653274095082324^-97.60372292227069 == 4.088393948750035e279
     @test 8.758520413376658e-5^70.55863059215994 == 5.052076767078296e-287
+
+    # issue #53881
+    @test 2.0 ^ typemin(Int) == 0.0
+    @test (-1.0) ^ typemin(Int) == 1.0
+    Z = Int64(2)
+    @test prevfloat(1.0) ^ (-Z^54) ≈ 7.38905609893065
+    @test prevfloat(1.0) ^ (-Z^62) ≈ 2.2844135865231613e222
+    @test prevfloat(1.0) ^ (-Z^63) == Inf
+    @test prevfloat(1.0) ^ (Z^62-1) * prevfloat(1.0) ^ (-Z^62+1) == 1.0
+    n, x = -1065564664, 0.9999997040311492
+    @test abs(x^n - Float64(big(x)^n)) / eps(x^n) == 0
 end
 
 # Test that sqrt behaves correctly and doesn't exhibit fp80 double rounding.