diff options
author | upstream source tree <ports@midipix.org> | 2015-03-15 20:14:05 -0400 |
---|---|---|
committer | upstream source tree <ports@midipix.org> | 2015-03-15 20:14:05 -0400 |
commit | 554fd8c5195424bdbcabf5de30fdc183aba391bd (patch) | |
tree | 976dc5ab7fddf506dadce60ae936f43f58787092 /libgo/go/math/pow.go | |
download | cbb-gcc-4.6.4-15d2061ac0796199866debe9ac87130894b0cdd3.tar.bz2 cbb-gcc-4.6.4-15d2061ac0796199866debe9ac87130894b0cdd3.tar.xz |
obtained gcc-4.6.4.tar.bz2 from upstream website;upstream
verified gcc-4.6.4.tar.bz2.sig;
imported gcc-4.6.4 source tree from verified upstream tarball.
downloading a git-generated archive based on the 'upstream' tag
should provide you with a source tree that is binary identical
to the one extracted from the above tarball.
if you have obtained the source via the command 'git clone',
however, do note that line-endings of files in your working
directory might differ from line-endings of the respective
files in the upstream repository.
Diffstat (limited to 'libgo/go/math/pow.go')
-rw-r--r-- | libgo/go/math/pow.go | 139 |
1 files changed, 139 insertions, 0 deletions
diff --git a/libgo/go/math/pow.go b/libgo/go/math/pow.go new file mode 100644 index 000000000..06b107401 --- /dev/null +++ b/libgo/go/math/pow.go @@ -0,0 +1,139 @@ +// 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 { + // TODO(rsc): Remove manual inlining of IsNaN, IsInf + // when compiler does it for us + 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 x != x || y != y: // 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 y > MaxFloat64 || y < -MaxFloat64: // IsInf(y, 0): + switch { + case x == -1: + return 1 + case (Fabs(x) < 1) == IsInf(y, 1): + return 0 + default: + return Inf(1) + } + case x > MaxFloat64 || x < -MaxFloat64: // 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) +} |