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.
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