math: eliminate overflow in Pow(x,y) for large y

The current implementation uses a shift and add
loop to compute the product of x's exponent xe and
the integer part of y (yi) for yi up to 1<<63.
Since xe is an 11-bit exponent, this product can be
up to 74-bits and overflow both 32 and 64-bit int.

This change checks whether the accumulated exponent
will fit in the 11-bit float exponent of the output
and breaks out of the loop early if overflow is detected.

The current handling of yi >= 1<<63 uses Exp(y * Log(x))
which incorrectly returns Nan for x<0.  In addition,
for y this large, Exp(y * Log(x)) can be enumerated
to only overflow except when x == -1 since the
boundary cases computed exactly:

Pow(NextAfter(1.0, Inf(1)), 1<<63)  == 2.72332... * 10^889
Pow(NextAfter(1.0, Inf(-1)), 1<<63) == 1.91624... * 10^-445

exceed the range of float64. So, the call can be
replaced with a simple case statement analgous to
y == Inf that correctly handles x < 0 as well.

Fixes #7394

Change-Id: I6f50dc951f3693697f9669697599860604323102
Reviewed-on: https://go-review.googlesource.com/48290
Reviewed-by: Robert Griesemer <gri@golang.org>

diff --git a/src/math/all_test.go b/src/math/all_test.go
index 4449228c1ea777a0238ad7ec4eca34972af54f2c..bdc4d228d59ae54af58cbd82517e3a3809f2c175 100644 (file)
--- a/src/math/all_test.go
+++ b/src/math/all_test.go
@@ -1586,6 +1586,17 @@ var vfpowSC = [][2]float64{
        {NaN(), 1},
        {NaN(), Pi},
        {NaN(), NaN()},
+
+       // Issue #7394 overflow checks
+       {2, float64(1 << 32)},
+       {2, -float64(1 << 32)},
+       {-2, float64(1<<32 + 1)},
+       {1 / 2, float64(1 << 45)},
+       {1 / 2, -float64(1 << 45)},
+       {Nextafter(1, 2), float64(1 << 63)},
+       {Nextafter(1, -2), float64(1 << 63)},
+       {Nextafter(-1, 2), float64(1 << 63)},
+       {Nextafter(-1, -2), float64(1 << 63)},
 }
 var powSC = []float64{
        0,               // pow(-Inf, -Pi)
@@ -1647,6 +1658,17 @@ var powSC = []float64{
        NaN(),           // pow(NaN, 1)
        NaN(),           // pow(NaN, +Pi)
        NaN(),           // pow(NaN, NaN)
+
+       // Issue #7394 overflow checks
+       Inf(1),  // pow(2, float64(1 << 32))
+       0,       // pow(2, -float64(1 << 32))
+       Inf(-1), // pow(-2, float64(1<<32 + 1))
+       0,       // pow(1/2, float64(1 << 45))
+       Inf(1),  // pow(1/2, -float64(1 << 45))
+       Inf(1),  // pow(Nextafter(1, 2), float64(1 << 63))
+       0,       // pow(Nextafter(1, -2), float64(1 << 63))
+       0,       // pow(Nextafter(-1, 2), float64(1 << 63))
+       Inf(1),  // pow(Nextafter(-1, -2), float64(1 << 63))
 }
 
 var vfpow10SC = []int{

diff --git a/src/math/pow.go b/src/math/pow.go
index b3bfadfb875a3dea456314591f64db28453a4fe2..daebf947284047f63d0a14fc7203dfab2af60f07 100644 (file)
--- a/src/math/pow.go
+++ b/src/math/pow.go
@@ -94,7 +94,16 @@ func pow(x, y float64) float64 {
                return NaN()
        }
        if yi >= 1<<63 {
-               return Exp(y * Log(x))
+               // yi is a large even int that will lead to overflow (or underflow to 0)
+               // for all x except -1 (x == 1 was handled earlier)
+               switch {
+               case x == -1:
+                       return 1
+               case (Abs(x) < 1) == (y > 0):
+                       return 0
+               default:
+                       return Inf(1)
+               }
        }
 
        // ans = a1 * 2**ae (= 1 for now).
@@ -116,6 +125,15 @@ func pow(x, y float64) float64 {
        // accumulate powers of two into exp.
        x1, xe := Frexp(x)
        for i := int64(yi); i != 0; i >>= 1 {
+               if xe < -1<<12 || 1<<12 < xe {
+                       // catch xe before it overflows the left shift below
+                       // Since i !=0 it has at least one bit still set, so ae will accumulate xe
+                       // on at least one more iteration, ae += xe is a lower bound on ae
+                       // the lower bound on ae exceeds the size of a float64 exp
+                       // so the final call to Ldexp will produce under/overflow (0/Inf)
+                       ae += xe
+                       break
+               }
                if i&1 == 1 {
                        a1 *= x1
                        ae += xe