1 // Copyright 2009 The Go Authors. All rights reserved.
2 // Use of this source code is governed by a BSD-style
3 // license that can be found in the LICENSE file.
7 func isOddInt(x
float64) bool {
9 return xf
== 0 && int64(xi
)&1 == 1
12 // Special cases taken from FreeBSD's /usr/src/lib/msun/src/e_pow.c
13 // updated by IEEE Std. 754-2008 "Section 9.2.1 Special values".
15 // Pow returns x**y, the base-x exponential of y.
17 // Special cases are (in order):
18 // Pow(x, ±0) = 1 for any x
19 // Pow(1, y) = 1 for any y
20 // Pow(x, 1) = x for any x
23 // Pow(±0, y) = ±Inf for y an odd integer < 0
24 // Pow(±0, -Inf) = +Inf
26 // Pow(±0, y) = +Inf for finite y < 0 and not an odd integer
27 // Pow(±0, y) = ±0 for y an odd integer > 0
28 // Pow(±0, y) = +0 for finite y > 0 and not an odd integer
30 // Pow(x, +Inf) = +Inf for |x| > 1
31 // Pow(x, -Inf) = +0 for |x| > 1
32 // Pow(x, +Inf) = +0 for |x| < 1
33 // Pow(x, -Inf) = +Inf for |x| < 1
34 // Pow(+Inf, y) = +Inf for y > 0
35 // Pow(+Inf, y) = +0 for y < 0
36 // Pow(-Inf, y) = Pow(-0, -y)
37 // Pow(x, y) = NaN for finite x < 0 and finite non-integer y
38 func Pow(x
, y
float64) float64 {
43 func libc_pow(float64, float64) float64
45 func pow(x
, y
float64) float64 {
47 case y
== 0 || x
== 1:
51 case IsNaN(x
) ||
IsNaN(y
):
57 return Copysign(Inf(1), x
)
70 case (Abs(x
) < 1) == IsInf(y
, 1):
77 return Pow(1/x
, -y
) // Pow(-0, -y)
91 yi
, yf
:= Modf(Abs(y
))
96 // yi is a large even int that will lead to overflow (or underflow to 0)
97 // for all x except -1 (x == 1 was handled earlier)
101 case (Abs(x
) < 1) == (y
> 0):
108 // ans = a1 * 2**ae (= 1 for now).
118 a1
= Exp(yf
* Log(x
))
122 // by multiplying in successive squarings
123 // of x according to bits of yi.
124 // accumulate powers of two into exp.
126 for i
:= int64(yi
); i
!= 0; i
>>= 1 {
127 if xe
< -1<<12 ||
1<<12 < xe
{
128 // catch xe before it overflows the left shift below
129 // Since i !=0 it has at least one bit still set, so ae will accumulate xe
130 // on at least one more iteration, ae += xe is a lower bound on ae
131 // the lower bound on ae exceeds the size of a float64 exp
132 // so the final call to Ldexp will produce under/overflow (0/Inf)
149 // if y < 0 { ans = 1 / ans }
150 // but in the opposite order