# 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  	xi, xf := Modf(x)
9  	return xf == 0 && int64(xi)&1 == 1
10  }
11
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".
14
15  // Pow returns x**y, the base-x exponential of y.
16  //
17  // Special cases are (in order):
18  //
19  //	Pow(x, ±0) = 1 for any x
20  //	Pow(1, y) = 1 for any y
21  //	Pow(x, 1) = x for any x
22  //	Pow(NaN, y) = NaN
23  //	Pow(x, NaN) = NaN
24  //	Pow(±0, y) = ±Inf for y an odd integer < 0
25  //	Pow(±0, -Inf) = +Inf
26  //	Pow(±0, +Inf) = +0
27  //	Pow(±0, y) = +Inf for finite y < 0 and not an odd integer
28  //	Pow(±0, y) = ±0 for y an odd integer > 0
29  //	Pow(±0, y) = +0 for finite y > 0 and not an odd integer
30  //	Pow(-1, ±Inf) = 1
31  //	Pow(x, +Inf) = +Inf for |x| > 1
32  //	Pow(x, -Inf) = +0 for |x| > 1
33  //	Pow(x, +Inf) = +0 for |x| < 1
34  //	Pow(x, -Inf) = +Inf for |x| < 1
35  //	Pow(+Inf, y) = +Inf for y > 0
36  //	Pow(+Inf, y) = +0 for y < 0
37  //	Pow(-Inf, y) = Pow(-0, -y)
38  //	Pow(x, y) = NaN for finite x < 0 and finite non-integer y
39  func Pow(x, y float64) float64 {
40  	if haveArchPow {
41  		return archPow(x, y)
42  	}
43  	return pow(x, y)
44  }
45
46  func pow(x, y float64) float64 {
47  	switch {
48  	case y == 0 || x == 1:
49  		return 1
50  	case y == 1:
51  		return x
52  	case IsNaN(x) || IsNaN(y):
53  		return NaN()
54  	case x == 0:
55  		switch {
56  		case y < 0:
57  			if isOddInt(y) {
58  				return Copysign(Inf(1), x)
59  			}
60  			return Inf(1)
61  		case y > 0:
62  			if isOddInt(y) {
63  				return x
64  			}
65  			return 0
66  		}
67  	case IsInf(y, 0):
68  		switch {
69  		case x == -1:
70  			return 1
71  		case (Abs(x) < 1) == IsInf(y, 1):
72  			return 0
73  		default:
74  			return Inf(1)
75  		}
76  	case IsInf(x, 0):
77  		if IsInf(x, -1) {
78  			return Pow(1/x, -y) // Pow(-0, -y)
79  		}
80  		switch {
81  		case y < 0:
82  			return 0
83  		case y > 0:
84  			return Inf(1)
85  		}
86  	case y == 0.5:
87  		return Sqrt(x)
88  	case y == -0.5:
89  		return 1 / Sqrt(x)
90  	}
91
92  	yi, yf := Modf(Abs(y))
93  	if yf != 0 && x < 0 {
94  		return NaN()
95  	}
96  	if yi >= 1<<63 {
97  		// yi is a large even int that will lead to overflow (or underflow to 0)
98  		// for all x except -1 (x == 1 was handled earlier)
99  		switch {
100  		case x == -1:
101  			return 1
102  		case (Abs(x) < 1) == (y > 0):
103  			return 0
104  		default:
105  			return Inf(1)
106  		}
107  	}
108
109  	// ans = a1 * 2**ae (= 1 for now).
110  	a1 := 1.0
111  	ae := 0
112
113  	// ans *= x**yf
114  	if yf != 0 {
115  		if yf > 0.5 {
116  			yf--
117  			yi++
118  		}
119  		a1 = Exp(yf * Log(x))
120  	}
121
122  	// ans *= x**yi
123  	// by multiplying in successive squarings
124  	// of x according to bits of yi.
125  	// accumulate powers of two into exp.
126  	x1, xe := Frexp(x)
127  	for i := int64(yi); i != 0; i >>= 1 {
128  		if xe < -1<<12 || 1<<12 < xe {
129  			// catch xe before it overflows the left shift below
130  			// Since i !=0 it has at least one bit still set, so ae will accumulate xe
131  			// on at least one more iteration, ae += xe is a lower bound on ae
132  			// the lower bound on ae exceeds the size of a float64 exp
133  			// so the final call to Ldexp will produce under/overflow (0/Inf)
134  			ae += xe
135  			break
136  		}
137  		if i&1 == 1 {
138  			a1 *= x1
139  			ae += xe
140  		}
141  		x1 *= x1
142  		xe <<= 1
143  		if x1 < .5 {
144  			x1 += x1
145  			xe--
146  		}
147  	}
148
149  	// ans = a1*2**ae
150  	// if y < 0 { ans = 1 / ans }
151  	// but in the opposite order
152  	if y < 0 {
153  		a1 = 1 / a1
154  		ae = -ae
155  	}
156  	return Ldexp(a1, ae)
157  }
158
```

View as plain text