Source file src/math/pow.go

```     1  // Copyright 2009 The Go Authors. All rights reserved.
2  // Use of this source code is governed by a BSD-style
4
5  package math
6
7  func isOddInt(x float64) bool {
8  	if Abs(x) >= (1 << 53) {
9  		// 1 << 53 is the largest exact integer in the float64 format.
10  		// Any number outside this range will be truncated before the decimal point and therefore will always be
11  		// an even integer.
12  		// Without this check and if x overflows int64 the int64(xi) conversion below may produce incorrect results
13  		// on some architectures (and does so on arm64). See issue #57465.
14  		return false
15  	}
16
17  	xi, xf := Modf(x)
18  	return xf == 0 && int64(xi)&1 == 1
19  }
20
21  // Special cases taken from FreeBSD's /usr/src/lib/msun/src/e_pow.c
22  // updated by IEEE Std. 754-2008 "Section 9.2.1 Special values".
23
24  // Pow returns x**y, the base-x exponential of y.
25  //
26  // Special cases are (in order):
27  //
28  //	Pow(x, ±0) = 1 for any x
29  //	Pow(1, y) = 1 for any y
30  //	Pow(x, 1) = x for any x
31  //	Pow(NaN, y) = NaN
32  //	Pow(x, NaN) = NaN
33  //	Pow(±0, y) = ±Inf for y an odd integer < 0
34  //	Pow(±0, -Inf) = +Inf
35  //	Pow(±0, +Inf) = +0
36  //	Pow(±0, y) = +Inf for finite y < 0 and not an odd integer
37  //	Pow(±0, y) = ±0 for y an odd integer > 0
38  //	Pow(±0, y) = +0 for finite y > 0 and not an odd integer
39  //	Pow(-1, ±Inf) = 1
40  //	Pow(x, +Inf) = +Inf for |x| > 1
41  //	Pow(x, -Inf) = +0 for |x| > 1
42  //	Pow(x, +Inf) = +0 for |x| < 1
43  //	Pow(x, -Inf) = +Inf for |x| < 1
44  //	Pow(+Inf, y) = +Inf for y > 0
45  //	Pow(+Inf, y) = +0 for y < 0
46  //	Pow(-Inf, y) = Pow(-0, -y)
47  //	Pow(x, y) = NaN for finite x < 0 and finite non-integer y
48  func Pow(x, y float64) float64 {
49  	if haveArchPow {
50  		return archPow(x, y)
51  	}
52  	return pow(x, y)
53  }
54
55  func pow(x, y float64) float64 {
56  	switch {
57  	case y == 0 || x == 1:
58  		return 1
59  	case y == 1:
60  		return x
61  	case IsNaN(x) || IsNaN(y):
62  		return NaN()
63  	case x == 0:
64  		switch {
65  		case y < 0:
66  			if Signbit(x) && isOddInt(y) {
67  				return Inf(-1)
68  			}
69  			return Inf(1)
70  		case y > 0:
71  			if Signbit(x) && isOddInt(y) {
72  				return x
73  			}
74  			return 0
75  		}
76  	case IsInf(y, 0):
77  		switch {
78  		case x == -1:
79  			return 1
80  		case (Abs(x) < 1) == IsInf(y, 1):
81  			return 0
82  		default:
83  			return Inf(1)
84  		}
85  	case IsInf(x, 0):
86  		if IsInf(x, -1) {
87  			return Pow(1/x, -y) // Pow(-0, -y)
88  		}
89  		switch {
90  		case y < 0:
91  			return 0
92  		case y > 0:
93  			return Inf(1)
94  		}
95  	case y == 0.5:
96  		return Sqrt(x)
97  	case y == -0.5:
98  		return 1 / Sqrt(x)
99  	}
100
101  	yi, yf := Modf(Abs(y))
102  	if yf != 0 && x < 0 {
103  		return NaN()
104  	}
105  	if yi >= 1<<63 {
106  		// yi is a large even int that will lead to overflow (or underflow to 0)
107  		// for all x except -1 (x == 1 was handled earlier)
108  		switch {
109  		case x == -1:
110  			return 1
111  		case (Abs(x) < 1) == (y > 0):
112  			return 0
113  		default:
114  			return Inf(1)
115  		}
116  	}
117
118  	// ans = a1 * 2**ae (= 1 for now).
119  	a1 := 1.0
120  	ae := 0
121
122  	// ans *= x**yf
123  	if yf != 0 {
124  		if yf > 0.5 {
125  			yf--
126  			yi++
127  		}
128  		a1 = Exp(yf * Log(x))
129  	}
130
131  	// ans *= x**yi
132  	// by multiplying in successive squarings
133  	// of x according to bits of yi.
134  	// accumulate powers of two into exp.
135  	x1, xe := Frexp(x)
136  	for i := int64(yi); i != 0; i >>= 1 {
137  		if xe < -1<<12 || 1<<12 < xe {
138  			// catch xe before it overflows the left shift below
139  			// Since i !=0 it has at least one bit still set, so ae will accumulate xe
140  			// on at least one more iteration, ae += xe is a lower bound on ae
141  			// the lower bound on ae exceeds the size of a float64 exp
142  			// so the final call to Ldexp will produce under/overflow (0/Inf)
143  			ae += xe
144  			break
145  		}
146  		if i&1 == 1 {
147  			a1 *= x1
148  			ae += xe
149  		}
150  		x1 *= x1
151  		xe <<= 1
152  		if x1 < .5 {
153  			x1 += x1
154  			xe--
155  		}
156  	}
157
158  	// ans = a1*2**ae
159  	// if y < 0 { ans = 1 / ans }
160  	// but in the opposite order
161  	if y < 0 {
162  		a1 = 1 / a1
163  		ae = -ae
164  	}
165  	return Ldexp(a1, ae)
166  }
167
```

View as plain text