blob: 336193bce12bf194eeb1152c620063911c7c85e6 [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 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) } 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 { // 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 flip { ans = 1 / ans } // but in the opposite order if flip { a1 = 1 / a1 ae = -ae } return Ldexp(a1, ae) }