138 lines
2.6 KiB
Go
138 lines
2.6 KiB
Go
// 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 {
|
|
switch {
|
|
case y == 0 || x == 1:
|
|
return 1
|
|
case y == 1:
|
|
return x
|
|
case y == 0.5:
|
|
return Sqrt(x)
|
|
case y == -0.5:
|
|
return 1 / Sqrt(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)
|
|
}
|
|
}
|
|
|
|
absy := y
|
|
flip := false
|
|
if absy < 0 {
|
|
absy = -absy
|
|
flip = true
|
|
}
|
|
yi, yf := Modf(absy)
|
|
if yf != 0 && x < 0 {
|
|
return NaN()
|
|
}
|
|
if yi >= 1<<63 {
|
|
return Exp(y * Log(x))
|
|
}
|
|
|
|
// 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 i&1 == 1 {
|
|
a1 *= x1
|
|
ae += xe
|
|
}
|
|
x1 *= x1
|
|
xe <<= 1
|
|
if x1 < .5 {
|
|
x1 += x1
|
|
xe--
|
|
}
|
|
}
|
|
|
|
// ans = a1*2**ae
|
|
// if flip { ans = 1 / ans }
|
|
// but in the opposite order
|
|
if flip {
|
|
a1 = 1 / a1
|
|
ae = -ae
|
|
}
|
|
return Ldexp(a1, ae)
|
|
}
|