blob: e45a044ae1cc02a2ae2d696579c6d1968b4c7b2c [file] [log] [blame]
// Copyright 2009 The Go Authors. All rights reserved.
// Use of this source code is governed by a BSD-style
// license that can be found in the LICENSE file.
package math
func isOddInt(x float64) bool {
xi, xf := Modf(x)
return xf == 0 && int64(xi)&1 == 1
}
// Special cases taken from FreeBSD's /usr/src/lib/msun/src/e_pow.c
// updated by IEEE Std. 754-2008 "Section 9.2.1 Special values".
// Pow returns x**y, the base-x exponential of y.
//
// Special cases are (in order):
// Pow(x, ±0) = 1 for any x
// Pow(1, y) = 1 for any y
// Pow(x, 1) = x for any x
// Pow(NaN, y) = NaN
// Pow(x, NaN) = NaN
// Pow(±0, y) = ±Inf for y an odd integer < 0
// Pow(±0, -Inf) = +Inf
// Pow(±0, +Inf) = +0
// Pow(±0, y) = +Inf for finite y < 0 and not an odd integer
// Pow(±0, y) = ±0 for y an odd integer > 0
// Pow(±0, y) = +0 for finite y > 0 and not an odd integer
// Pow(-1, ±Inf) = 1
// Pow(x, +Inf) = +Inf for |x| > 1
// Pow(x, -Inf) = +0 for |x| > 1
// Pow(x, +Inf) = +0 for |x| < 1
// Pow(x, -Inf) = +Inf for |x| < 1
// Pow(+Inf, y) = +Inf for y > 0
// Pow(+Inf, y) = +0 for y < 0
// Pow(-Inf, y) = Pow(-0, -y)
// Pow(x, y) = NaN for finite x < 0 and finite non-integer y
func Pow(x, y float64) float64 {
if haveArchPow {
return archPow(x, y)
}
return pow(x, y)
}
func pow(x, y float64) float64 {
switch {
case y == 0 || x == 1:
return 1
case y == 1:
return x
case IsNaN(x) || IsNaN(y):
return NaN()
case x == 0:
switch {
case y < 0:
if isOddInt(y) {
return Copysign(Inf(1), x)
}
return Inf(1)
case y > 0:
if isOddInt(y) {
return x
}
return 0
}
case IsInf(y, 0):
switch {
case x == -1:
return 1
case (Abs(x) < 1) == IsInf(y, 1):
return 0
default:
return Inf(1)
}
case IsInf(x, 0):
if IsInf(x, -1) {
return Pow(1/x, -y) // Pow(-0, -y)
}
switch {
case y < 0:
return 0
case y > 0:
return Inf(1)
}
case y == 0.5:
return Sqrt(x)
case y == -0.5:
return 1 / Sqrt(x)
}
yi, yf := Modf(Abs(y))
if yf != 0 && x < 0 {
return NaN()
}
if yi >= 1<<63 {
// 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).
a1 := 1.0
ae := 0
// ans *= x**yf
if yf != 0 {
if yf > 0.5 {
yf--
yi++
}
a1 = Exp(yf * Log(x))
}
// ans *= x**yi
// by multiplying in successive squarings
// of x according to bits of yi.
// 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
}
x1 *= x1
xe <<= 1
if x1 < .5 {
x1 += x1
xe--
}
}
// ans = a1*2**ae
// if y < 0 { ans = 1 / ans }
// but in the opposite order
if y < 0 {
a1 = 1 / a1
ae = -ae
}
return Ldexp(a1, ae)
}