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/48290Reviewed-by: Robert Griesemer <gri@golang.org>

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{

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