From 979027bf0be8f7900b528d394a3170bb198a281c Mon Sep 17 00:00:00 2001
From: KlausC <klausC@users.noreply.github.com>
Date: Thu, 28 Mar 2024 19:42:33 +0100
Subject: [PATCH 01/14] fix `Math.pow_body` for huge powers

---
 base/math.jl | 27 ++++++++++++++++++++++-----
 test/math.jl |  7 +++++++
 2 files changed, 29 insertions(+), 5 deletions(-)

diff --git a/base/math.jl b/base/math.jl
index c3c15505c5f8f..99561cf41ffbe 100644
--- a/base/math.jl
+++ b/base/math.jl
@@ -1316,17 +1316,34 @@ end
 
 @assume_effects :terminates_locally @noinline function pow_body(x::Float64, n::Integer)
     y = 1.0
-    xnlo = ynlo = 0.0
+    xnlo = ynlo = err = 0.0
+    if iszero(n)
+        return y
+    elseif n == 1
+        return x
+    elseif n == -1
+        return inv(x)
+    end
+    x0 = x
     n == 3 && return x*x*x # keep compatibility with literal_pow
+    if isodd(n)
+        y = n < 0 ? inv(x) : x
+    end
+    x, xnlo = two_mul(x, x)
+    m = n ÷ 2
     if n < 0
         rx = inv(x)
-        n==-2 && return rx*rx #keep compatibility with literal_pow
-        isfinite(x) && (xnlo = -fma(x, rx, -1.) * rx)
+        n == -2 && return rx #keep compatibility with literal_pow
+        !isfinite(rx) && return isodd(n) ? copysign(rx, x0) : rx
+        if isfinite(x)
+            xnlo = -(xnlo * rx + fma(x, rx, -1.0)) * rx
+        end
         x = rx
-        n = -n
+        m = -m
     end
+    n = m
     while n > 1
-        if n&1 > 0
+        if isodd(n)
             err = muladd(y, xnlo, x*ynlo)
             y, ynlo = two_mul(x,y)
             ynlo += err
diff --git a/test/math.jl b/test/math.jl
index f551bb3a5d4b7..e369d81a08d22 100644
--- a/test/math.jl
+++ b/test/math.jl
@@ -1459,6 +1459,13 @@ end
     # 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
+    @test prevfloat(1.0) ^ (-2^54) ≈ 7.38905609893065
+    @test prevfloat(1.0) ^ (-2^62) ≈ 2.2844135865231613e222
+    @test prevfloat(1.0) ^ (-2^63) == Inf
 end
 
 # Test that sqrt behaves correctly and doesn't exhibit fp80 double rounding.

From ed67f6a93e1eca1c1b74a0a769d22adb7fa6bb37 Mon Sep 17 00:00:00 2001
From: KlausC <klausC@users.noreply.github.com>
Date: Thu, 28 Mar 2024 21:17:39 +0100
Subject: [PATCH 02/14] add compensation in `pow_body`

---
 base/math.jl | 7 ++++++-
 1 file changed, 6 insertions(+), 1 deletion(-)

diff --git a/base/math.jl b/base/math.jl
index 99561cf41ffbe..7d3da1f68c267 100644
--- a/base/math.jl
+++ b/base/math.jl
@@ -1327,7 +1327,12 @@ end
     x0 = x
     n == 3 && return x*x*x # keep compatibility with literal_pow
     if isodd(n)
-        y = n < 0 ? inv(x) : x
+        if n < 0
+            y = inv(x)
+            ynlo = -fma(x, y, -1.0) * y
+        else
+            y = x
+        end
     end
     x, xnlo = two_mul(x, x)
     m = n ÷ 2

From da81b2b933b6bd5cf64397a3e95549c832d184f7 Mon Sep 17 00:00:00 2001
From: KlausC <klausC@users.noreply.github.com>
Date: Thu, 28 Mar 2024 21:31:14 +0100
Subject: [PATCH 03/14] cosmetics

---
 base/math.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/base/math.jl b/base/math.jl
index 7d3da1f68c267..fe18261c4eddb 100644
--- a/base/math.jl
+++ b/base/math.jl
@@ -1329,7 +1329,7 @@ end
     if isodd(n)
         if n < 0
             y = inv(x)
-            ynlo = -fma(x, y, -1.0) * y
+            ynlo = fma(-x, y, 1.0) * y
         else
             y = x
         end
@@ -1341,7 +1341,7 @@ end
         n == -2 && return rx #keep compatibility with literal_pow
         !isfinite(rx) && return isodd(n) ? copysign(rx, x0) : rx
         if isfinite(x)
-            xnlo = -(xnlo * rx + fma(x, rx, -1.0)) * rx
+            xnlo = (fma(-x, rx, 1.0) - xnlo * rx) * rx
         end
         x = rx
         m = -m

From df4292af7b0a8467e2df4f240f3fa36f0f6ed315 Mon Sep 17 00:00:00 2001
From: KlausC <klausC@users.noreply.github.com>
Date: Thu, 28 Mar 2024 21:40:20 +0100
Subject: [PATCH 04/14] another call to `fma`

---
 base/math.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/base/math.jl b/base/math.jl
index fe18261c4eddb..6225771c968cb 100644
--- a/base/math.jl
+++ b/base/math.jl
@@ -1341,7 +1341,7 @@ end
         n == -2 && return rx #keep compatibility with literal_pow
         !isfinite(rx) && return isodd(n) ? copysign(rx, x0) : rx
         if isfinite(x)
-            xnlo = (fma(-x, rx, 1.0) - xnlo * rx) * rx
+            xnlo = fma(-xnlo, rx, (fma(-x, rx, 1.0))) * rx
         end
         x = rx
         m = -m

From cd5dd728e3cf15e95088761a3e489a8df3cd9f1b Mon Sep 17 00:00:00 2001
From: KlausC <klausC@users.noreply.github.com>
Date: Fri, 29 Mar 2024 05:14:26 +0100
Subject: [PATCH 05/14] kahans sum - fix power -2

---
 base/math.jl | 13 +++++++++----
 1 file changed, 9 insertions(+), 4 deletions(-)

diff --git a/base/math.jl b/base/math.jl
index 6225771c968cb..84e5f78ab7314 100644
--- a/base/math.jl
+++ b/base/math.jl
@@ -1316,7 +1316,7 @@ end
 
 @assume_effects :terminates_locally @noinline function pow_body(x::Float64, n::Integer)
     y = 1.0
-    xnlo = ynlo = err = 0.0
+    xnlo = ynlo = 0.0
     if iszero(n)
         return y
     elseif n == 1
@@ -1337,8 +1337,11 @@ end
     x, xnlo = two_mul(x, x)
     m = n ÷ 2
     if n < 0
+        if n == -2
+            rx = inv(x0)
+            return rx * rx #keep compatibility with literal_pow
+        end
         rx = inv(x)
-        n == -2 && return rx #keep compatibility with literal_pow
         !isfinite(rx) && return isodd(n) ? copysign(rx, x0) : rx
         if isfinite(x)
             xnlo = fma(-xnlo, rx, (fma(-x, rx, 1.0))) * rx
@@ -1353,9 +1356,11 @@ end
             y, ynlo = two_mul(x,y)
             ynlo += err
         end
-        err = x*2*xnlo
-        x, xnlo = two_mul(x, x)
+        err = (xnlo + xnlo) * x
+        xx, xnlo = two_mul(x, x)
         xnlo += err
+        x = xx + xnlo
+        xnlo -= x - xx
         n >>>= 1
     end
     err = muladd(y, xnlo, x*ynlo)

From fc0b9f5784624dd5b3481f2a230eb8e4fb78ce7f Mon Sep 17 00:00:00 2001
From: KlausC <klausC@users.noreply.github.com>
Date: Fri, 29 Mar 2024 05:38:07 +0100
Subject: [PATCH 06/14] fix infinite case

---
 base/math.jl | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/base/math.jl b/base/math.jl
index 84e5f78ab7314..2984df794e80f 100644
--- a/base/math.jl
+++ b/base/math.jl
@@ -1324,7 +1324,7 @@ end
     elseif n == -1
         return inv(x)
     end
-    x0 = x
+    xx = x
     n == 3 && return x*x*x # keep compatibility with literal_pow
     if isodd(n)
         if n < 0
@@ -1338,11 +1338,11 @@ end
     m = n ÷ 2
     if n < 0
         if n == -2
-            rx = inv(x0)
+            rx = inv(xx)
             return rx * rx #keep compatibility with literal_pow
         end
         rx = inv(x)
-        !isfinite(rx) && return isodd(n) ? copysign(rx, x0) : rx
+        !isfinite(rx) && return isodd(n) ? copysign(rx, xx) : rx
         if isfinite(x)
             xnlo = fma(-xnlo, rx, (fma(-x, rx, 1.0))) * rx
         end
@@ -1364,7 +1364,7 @@ end
         n >>>= 1
     end
     err = muladd(y, xnlo, x*ynlo)
-    return ifelse(isfinite(x) & isfinite(err), muladd(x, y, err), x*y)
+    return ifelse(isfinite(x) & isfinite(err), muladd(x, y, err), xx * y)
 end
 
 function ^(x::Float32, n::Integer)

From 8cc2f9bc79c59da1b86d98bf5357dbb3faf23d23 Mon Sep 17 00:00:00 2001
From: KlausC <klausC@users.noreply.github.com>
Date: Fri, 29 Mar 2024 06:08:36 +0100
Subject: [PATCH 07/14] break on infinite

---
 base/math.jl | 5 +++--
 1 file changed, 3 insertions(+), 2 deletions(-)

diff --git a/base/math.jl b/base/math.jl
index 2984df794e80f..10938e01eae89 100644
--- a/base/math.jl
+++ b/base/math.jl
@@ -1337,9 +1337,9 @@ end
     x, xnlo = two_mul(x, x)
     m = n ÷ 2
     if n < 0
+        xx = inv(xx)
         if n == -2
-            rx = inv(xx)
-            return rx * rx #keep compatibility with literal_pow
+            return xx * xx #keep compatibility with literal_pow
         end
         rx = inv(x)
         !isfinite(rx) && return isodd(n) ? copysign(rx, xx) : rx
@@ -1361,6 +1361,7 @@ end
         xnlo += err
         x = xx + xnlo
         xnlo -= x - xx
+        isfinite(x) || break
         n >>>= 1
     end
     err = muladd(y, xnlo, x*ynlo)

From 5799dd90f976b7a215bcd7c2660aef53a8f42d84 Mon Sep 17 00:00:00 2001
From: KlausC <klausC@users.noreply.github.com>
Date: Fri, 29 Mar 2024 11:19:19 +0100
Subject: [PATCH 08/14] tests for x86 - use `muladd` for `fma`

---
 base/math.jl | 6 ++++--
 test/math.jl | 8 +++++---
 2 files changed, 9 insertions(+), 5 deletions(-)

diff --git a/base/math.jl b/base/math.jl
index 10938e01eae89..2ff298fb238f6 100644
--- a/base/math.jl
+++ b/base/math.jl
@@ -1314,6 +1314,8 @@ 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
@@ -1329,7 +1331,7 @@ end
     if isodd(n)
         if n < 0
             y = inv(x)
-            ynlo = fma(-x, y, 1.0) * y
+            ynlo = muladd(-x, y, 1.0) * y
         else
             y = x
         end
@@ -1344,7 +1346,7 @@ end
         rx = inv(x)
         !isfinite(rx) && return isodd(n) ? copysign(rx, xx) : rx
         if isfinite(x)
-            xnlo = fma(-xnlo, rx, (fma(-x, rx, 1.0))) * rx
+            xnlo = muladd(-xnlo, rx, (muladd(-x, rx, 1.0))) * rx
         end
         x = rx
         m = -m
diff --git a/test/math.jl b/test/math.jl
index e369d81a08d22..a800b0fcb55f3 100644
--- a/test/math.jl
+++ b/test/math.jl
@@ -1463,9 +1463,11 @@ end
     # issue #53881
     @test 2.0 ^ typemin(Int) == 0.0
     @test (-1.0) ^ typemin(Int) == 1.0
-    @test prevfloat(1.0) ^ (-2^54) ≈ 7.38905609893065
-    @test prevfloat(1.0) ^ (-2^62) ≈ 2.2844135865231613e222
-    @test prevfloat(1.0) ^ (-2^63) == Inf
+    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) * prevfloat(1.0) ^ (-Z^62) == 1.0
 end
 
 # Test that sqrt behaves correctly and doesn't exhibit fp80 double rounding.

From 7e9a370fc542fdc575ca2d5762920f74a1879637 Mon Sep 17 00:00:00 2001
From: KlausC <klausC@users.noreply.github.com>
Date: Fri, 29 Mar 2024 14:22:55 +0100
Subject: [PATCH 09/14] revert: use fma for muladd

---
 base/math.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/base/math.jl b/base/math.jl
index 2ff298fb238f6..ae5d68a8bd540 100644
--- a/base/math.jl
+++ b/base/math.jl
@@ -1331,7 +1331,7 @@ end
     if isodd(n)
         if n < 0
             y = inv(x)
-            ynlo = muladd(-x, y, 1.0) * y
+            ynlo = fma(-x, y, 1.0) * y
         else
             y = x
         end
@@ -1346,7 +1346,7 @@ end
         rx = inv(x)
         !isfinite(rx) && return isodd(n) ? copysign(rx, xx) : rx
         if isfinite(x)
-            xnlo = muladd(-xnlo, rx, (muladd(-x, rx, 1.0))) * rx
+            xnlo = fma(-xnlo, rx, (fma(-x, rx, 1.0))) * rx
         end
         x = rx
         m = -m

From 94749c7f838ed7c98fc6c59ff39659e2fb65dba2 Mon Sep 17 00:00:00 2001
From: KlausC <klausC@users.noreply.github.com>
Date: Sat, 30 Mar 2024 23:31:03 +0100
Subject: [PATCH 10/14] efficiency

---
 base/math.jl | 78 +++++++++++++++++++++++++---------------------------
 1 file changed, 38 insertions(+), 40 deletions(-)

diff --git a/base/math.jl b/base/math.jl
index ae5d68a8bd540..7d051ae822c09 100644
--- a/base/math.jl
+++ b/base/math.jl
@@ -1318,58 +1318,56 @@ end
 # calculations are effectively in double precision using fma
 @assume_effects :terminates_locally @noinline function pow_body(x::Float64, n::Integer)
     y = 1.0
+    n == 0 && return y
+    isnan(x) && return x
     xnlo = ynlo = 0.0
-    if iszero(n)
-        return y
-    elseif n == 1
-        return x
-    elseif n == -1
-        return inv(x)
-    end
-    xx = x
     n == 3 && return x*x*x # keep compatibility with literal_pow
-    if isodd(n)
-        if n < 0
-            y = inv(x)
-            ynlo = fma(-x, y, 1.0) * y
-        else
-            y = x
-        end
-    end
-    x, xnlo = two_mul(x, x)
-    m = n ÷ 2
+    negate = false
+    invert = n < 0
     if n < 0
-        xx = inv(xx)
         if n == -2
-            return xx * xx #keep compatibility with literal_pow
+            rx = inv(x)
+            return rx * rx #keep compatibility with literal_pow
         end
-        rx = inv(x)
-        !isfinite(rx) && return isodd(n) ? copysign(rx, xx) : rx
-        if isfinite(x)
-            xnlo = fma(-xnlo, rx, (fma(-x, rx, 1.0))) * rx
+        negate = n & 1 > 0
+        n = -n
+    end
+    y, ynlo = pow_loop(x, xnlo, y, ynlo, n)
+    if isfinite(y) && !iszero(y)
+        if invert
+            x = inv(y)
+            err = fma(x, y, -1.0)
+            y = (- ynlo * x + err) * x + x
         end
-        x = rx
-        m = -m
+    else
+        y = invert == iszero(y) ? Inf : 0.0
     end
-    n = m
-    while n > 1
-        if isodd(n)
-            err = muladd(y, xnlo, x*ynlo)
-            y, ynlo = two_mul(x,y)
-            ynlo += err
+    return copysign(y, negate)
+end
+
+function pow_loop(x, xnlo, y, ynlo, n)
+    for i = 1:sizeof(n)*8
+        if n & 1 > 0
+            err = y * xnlo + x * ynlo
+            xx = y
+            y = x * xx
+            ynlo = fma(x, xx, -y) + err
         end
-        err = (xnlo + xnlo) * x
-        xx, xnlo = two_mul(x, x)
-        xnlo += err
-        x = xx + xnlo
-        xnlo -= x - xx
-        isfinite(x) || break
         n >>>= 1
+        n == 0 && break
+        xx = x * x
+        xnlo = fma(x, x, -xx) + x * 2 * xnlo
+        if i == 30
+            x = xx + xnlo
+            xnlo -= x - xx
+        else
+            x = xx
+        end
     end
-    err = muladd(y, xnlo, x*ynlo)
-    return ifelse(isfinite(x) & isfinite(err), muladd(x, y, err), xx * y)
+    y, ynlo
 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

From 913d13fd3bfeeb0e83274c46b15362f72ecdff96 Mon Sep 17 00:00:00 2001
From: KlausC <klausC@users.noreply.github.com>
Date: Mon, 1 Apr 2024 12:44:08 +0200
Subject: [PATCH 11/14] accuacy and performance tweaked `pow_body`

---
 base/math.jl | 22 ++++++++++------------
 1 file changed, 10 insertions(+), 12 deletions(-)

diff --git a/base/math.jl b/base/math.jl
index 7d051ae822c09..2298a82a193db 100644
--- a/base/math.jl
+++ b/base/math.jl
@@ -1318,7 +1318,6 @@ end
 # calculations are effectively in double precision using fma
 @assume_effects :terminates_locally @noinline function pow_body(x::Float64, n::Integer)
     y = 1.0
-    n == 0 && return y
     isnan(x) && return x
     xnlo = ynlo = 0.0
     n == 3 && return x*x*x # keep compatibility with literal_pow
@@ -1332,12 +1331,16 @@ end
         negate = n & 1 > 0
         n = -n
     end
-    y, ynlo = pow_loop(x, xnlo, y, ynlo, n)
+    while n != 0
+        xx, xnlo, yy, ynlo, n = pow_loop(x, xnlo, y, ynlo, n, 16)
+        x = xx + xnlo; xnlo -= x - xx
+        y = yy + ynlo; ynlo -= y - yy
+    end
     if isfinite(y) && !iszero(y)
         if invert
             x = inv(y)
             err = fma(x, y, -1.0)
-            y = (- ynlo * x + err) * x + x
+            y = (-ynlo * x + err) * x + x
         end
     else
         y = invert == iszero(y) ? Inf : 0.0
@@ -1345,8 +1348,8 @@ end
     return copysign(y, negate)
 end
 
-function pow_loop(x, xnlo, y, ynlo, n)
-    for i = 1:sizeof(n)*8
+@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
@@ -1357,14 +1360,9 @@ function pow_loop(x, xnlo, y, ynlo, n)
         n == 0 && break
         xx = x * x
         xnlo = fma(x, x, -xx) + x * 2 * xnlo
-        if i == 30
-            x = xx + xnlo
-            xnlo -= x - xx
-        else
-            x = xx
-        end
+        x = xx
     end
-    y, ynlo
+    x, xnlo, y, ynlo, n
 end
 
 

From f43cdde49590db91762336ac621460003c817d9d Mon Sep 17 00:00:00 2001
From: KlausC <klausC@users.noreply.github.com>
Date: Mon, 1 Apr 2024 18:51:51 +0200
Subject: [PATCH 12/14] fixed tests and sign/infinity cases

---
 base/math.jl | 23 +++++++++++++++--------
 test/math.jl |  2 +-
 2 files changed, 16 insertions(+), 9 deletions(-)

diff --git a/base/math.jl b/base/math.jl
index 2298a82a193db..0cfcceff6a808 100644
--- a/base/math.jl
+++ b/base/math.jl
@@ -1319,17 +1319,26 @@ end
 @assume_effects :terminates_locally @noinline function pow_body(x::Float64, n::Integer)
     y = 1.0
     isnan(x) && return x
-    xnlo = ynlo = 0.0
     n == 3 && return x*x*x # keep compatibility with literal_pow
-    negate = false
-    invert = n < 0
+    xnlo = ynlo = 0.0
+    negate = signbit(x) && n & 1 > 0
+    xa = abs(x)
+    toinf = xa > 1
+    invert = false
     if n < 0
         if n == -2
             rx = inv(x)
             return rx * rx #keep compatibility with literal_pow
         end
-        negate = n & 1 > 0
+        if toinf
+            xx = inv(x)
+            xnlo = fma(-xx, x, 1.0) * xx
+            x = xx
+        else
+            invert = true
+        end
         n = -n
+        toinf = !toinf
     end
     while n != 0
         xx, xnlo, yy, ynlo, n = pow_loop(x, xnlo, y, ynlo, n, 16)
@@ -1342,13 +1351,11 @@ end
             err = fma(x, y, -1.0)
             y = (-ynlo * x + err) * x + x
         end
-    else
-        y = invert == iszero(y) ? Inf : 0.0
     end
-    return copysign(y, negate)
+    return isfinite(y) && !iszero(y) ? y : flipsign(toinf ? Inf : 0.0, -negate)
 end
 
-@inline function pow_loop(x, xnlo, y, ynlo, n, maxi)
+@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
diff --git a/test/math.jl b/test/math.jl
index a800b0fcb55f3..57f6d3e4fd86c 100644
--- a/test/math.jl
+++ b/test/math.jl
@@ -1454,7 +1454,7 @@ 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) ^ -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

From 788556c219e514dc1c85592c296dda2f8cb36ba8 Mon Sep 17 00:00:00 2001
From: KlausC <klausC@users.noreply.github.com>
Date: Mon, 1 Apr 2024 21:52:31 +0200
Subject: [PATCH 13/14] Test case for Int32

---
 test/math.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/math.jl b/test/math.jl
index 57f6d3e4fd86c..862aeeafd5dde 100644
--- a/test/math.jl
+++ b/test/math.jl
@@ -1454,7 +1454,7 @@ end
         @test nextfloat(T(-1))^floatmax(T) === T(0.0)
     end
     # test for large negative exponent where error compensation matters
-    @test prevfloat(1.0) ^ -2^62 == 2.2844135865398217e222
+    @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

From 31e810aa37025972290fc36811445fa858de48f0 Mon Sep 17 00:00:00 2001
From: KlausC <klausC@users.noreply.github.com>
Date: Wed, 3 Apr 2024 17:27:11 +0200
Subject: [PATCH 14/14] polish the algorithm - 100% <= 1 ULPs

---
 base/math.jl | 32 ++++++++++++++++++--------------
 test/math.jl |  4 +++-
 2 files changed, 21 insertions(+), 15 deletions(-)

diff --git a/base/math.jl b/base/math.jl
index 0cfcceff6a808..197f50ebd5946 100644
--- a/base/math.jl
+++ b/base/math.jl
@@ -1320,39 +1320,43 @@ end
     y = 1.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
+    toinf = xa >= 1
     invert = false
     if n < 0
         if n == -2
             rx = inv(x)
             return rx * rx #keep compatibility with literal_pow
         end
-        if toinf
-            xx = inv(x)
-            xnlo = fma(-xx, x, 1.0) * xx
-            x = xx
+        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
+    MAXI = 16
+    maxi = invert ? 1 : MAXI
     while n != 0
-        xx, xnlo, yy, ynlo, n = pow_loop(x, xnlo, y, ynlo, n, 16)
-        x = xx + xnlo; xnlo -= x - xx
-        y = yy + ynlo; ynlo -= y - yy
-    end
-    if isfinite(y) && !iszero(y)
+        x, xnlo = normal(x, xnlo)
+        x, xnlo, y, ynlo, n = pow_loop(x, xnlo, y, ynlo, n, maxi)
+
         if invert
-            x = inv(y)
-            err = fma(x, y, -1.0)
-            y = (-ynlo * x + err) * x + x
+            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(toinf ? Inf : 0.0, -negate)
+    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)
diff --git a/test/math.jl b/test/math.jl
index 862aeeafd5dde..179149e321996 100644
--- a/test/math.jl
+++ b/test/math.jl
@@ -1467,7 +1467,9 @@ end
     @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) * prevfloat(1.0) ^ (-Z^62) == 1.0
+    @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.