1// Ported from go, which is licensed under a BSD-3 license.
2// https://golang.org/LICENSE
3//
4// https://golang.org/src/math/pow.go
5
6const std = @import("../std.zig");
7const math = std.math;
8const expect = std.testing.expect;
9
10/// Returns x raised to the power of y (x^y).
11///
12/// Special Cases:
13/// - pow(x, +-0) = 1 for any x
14/// - pow(1, y) = 1 for any y
15/// - pow(x, 1) = x for any x
16/// - pow(nan, y) = nan for any y != 0
17/// - pow(x, nan) = nan for any x != 1
18/// - pow(+-0, y) = +-inf for y an odd integer < 0
19/// - pow(+-0, -inf) = +inf
20/// - pow(+-0, +inf) = +0
21/// - pow(+-0, y) = +inf for finite y < 0 and not an odd integer
22/// - pow(+-0, y) = +-0 for finite y an odd integer > 0
23/// - pow(+-0, y) = +0 for finite y > 0 and not an odd integer
24/// - pow(-1, +-inf) = 1
25/// - pow(x, +inf) = +inf for |x| > 1
26/// - pow(x, -inf) = +0 for |x| > 1
27/// - pow(x, +inf) = +0 for |x| < 1
28/// - pow(x, -inf) = +inf for |x| < 1
29/// - pow(+inf, y) = +inf for y > 0
30/// - pow(+inf, y) = +0 for y < 0
31/// - pow(-inf, y) = pow(-0, -y)
32/// - pow(x, y) = nan for finite x < 0 and finite non-integer y
33pub fn pow(comptime T: type, x: T, y: T) T {
34 if (@typeInfo(T) == .int) {
35 return math.powi(T, x, y) catch unreachable;
36 }
37
38 if (T != f32 and T != f64) {
39 @compileError("pow not implemented for " ++ @typeName(T));
40 }
41
42 // pow(x, +-0) = 1 for any x
43 // pow(1, y) = 1 for any y
44 if (y == 0 or x == 1) {
45 return 1;
46 }
47
48 // pow(nan, y) = nan for any y != 0
49 // pow(x, nan) = nan for any x != 1
50 if (math.isNan(x) or math.isNan(y)) {
51 @branchHint(.unlikely);
52 return math.nan(T);
53 }
54
55 // pow(x, 1) = x for any x
56 if (y == 1) {
57 return x;
58 }
59
60 if (x == 0) {
61 if (y < 0) {
62 // pow(+-0, y) = +-inf for y an odd integer
63 if (isOddInteger(y)) {
64 return math.copysign(math.inf(T), x);
65 }
66 // pow(+-0, y) = +inf for y an even integer
67 else {
68 return math.inf(T);
69 }
70 } else {
71 if (isOddInteger(y)) {
72 return x;
73 } else {
74 return 0;
75 }
76 }
77 }
78
79 if (math.isInf(y)) {
80 // pow(-1, inf) = 1
81 if (x == -1) {
82 return 1.0;
83 }
84 // pow(x, +inf) = +0 for |x| < 1
85 // pow(x, -inf) = +0 for |x| > 1
86 else if ((@abs(x) < 1) == math.isPositiveInf(y)) {
87 return 0;
88 }
89 // pow(x, -inf) = +inf for |x| < 1
90 // pow(x, +inf) = +inf for |x| > 1
91 else {
92 return math.inf(T);
93 }
94 }
95
96 if (math.isInf(x)) {
97 if (math.isNegativeInf(x)) {
98 return pow(T, 1 / x, -y);
99 }
100 // pow(+inf, y) = +0 for y < 0
101 else if (y < 0) {
102 return 0;
103 }
104 // pow(+inf, y) = +0 for y > 0
105 else if (y > 0) {
106 return math.inf(T);
107 }
108 }
109
110 // special case sqrt
111 if (y == 0.5) {
112 return @sqrt(x);
113 }
114
115 if (y == -0.5) {
116 return 1 / @sqrt(x);
117 }
118
119 const r1 = math.modf(@abs(y));
120 var yi = r1.ipart;
121 var yf = r1.fpart;
122
123 if (yf != 0 and x < 0) {
124 return math.nan(T);
125 }
126 if (yi >= 1 << (@typeInfo(T).float.bits - 1)) {
127 // yi is a large even int, so the result is always positive
128 // and the sign of x doesn't matter
129 return @exp(y * @log(@abs(x)));
130 }
131
132 // a = a1 * 2^ae
133 var a1: T = 1.0;
134 var ae: i32 = 0;
135
136 // a *= x^yf
137 if (yf != 0) {
138 if (yf > 0.5) {
139 yf -= 1;
140 yi += 1;
141 }
142 a1 = @exp(yf * @log(x));
143 }
144
145 // a *= x^yi
146 const r2 = math.frexp(x);
147 var xe = r2.exponent;
148 var x1 = r2.significand;
149
150 var i = @as(@Int(.signed, @typeInfo(T).float.bits), @intFromFloat(yi));
151 while (i != 0) : (i >>= 1) {
152 const overflow_shift = math.floatExponentBits(T) + 1;
153 if (xe < -(1 << overflow_shift) or (1 << overflow_shift) < xe) {
154 // catch xe before it overflows the left shift below
155 // Since i != 0 it has at least one bit still set, so ae will accumulate xe
156 // on at least one more iteration, ae += xe is a lower bound on ae
157 // the lower bound on ae exceeds the size of a float exp
158 // so the final call to Ldexp will produce under/overflow (0/Inf)
159 ae += xe;
160 break;
161 }
162 if (i & 1 == 1) {
163 a1 *= x1;
164 ae += xe;
165 }
166 x1 *= x1;
167 xe <<= 1;
168 if (x1 < 0.5) {
169 x1 += x1;
170 xe -= 1;
171 }
172 }
173
174 // a *= a1 * 2^ae
175 if (y < 0) {
176 a1 = 1 / a1;
177 ae = -ae;
178 }
179
180 return math.scalbn(a1, ae);
181}
182
183fn isOddInteger(x: f64) bool {
184 if (@abs(x) >= 1 << 53) {
185 // From https://golang.org/src/math/pow.go
186 // 1 << 53 is the largest exact integer in the float64 format.
187 // Any number outside this range will be truncated before the decimal point and therefore will always be
188 // an even integer.
189 // Without this check and if x overflows i64 the @intFromFloat(r.ipart) conversion below will panic
190 return false;
191 }
192 const r = math.modf(x);
193 return r.fpart == 0.0 and @as(i64, @intFromFloat(r.ipart)) & 1 == 1;
194}
195
196test isOddInteger {
197 try expect(isOddInteger(@floatFromInt(math.maxInt(i64) * 2)) == false);
198 try expect(isOddInteger(@floatFromInt(math.maxInt(i64) * 2 + 1)) == false);
199 try expect(isOddInteger(1 << 53) == false);
200 try expect(isOddInteger(12.0) == false);
201 try expect(isOddInteger(15.0) == true);
202}
203
204test pow {
205 const epsilon = 0.000001;
206
207 try expect(math.approxEqAbs(f32, pow(f32, 0.0, 3.3), 0.0, epsilon));
208 try expect(math.approxEqAbs(f32, pow(f32, 0.8923, 3.3), 0.686572, epsilon));
209 try expect(math.approxEqAbs(f32, pow(f32, 0.2, 3.3), 0.004936, epsilon));
210 try expect(math.approxEqAbs(f32, pow(f32, 1.5, 3.3), 3.811546, epsilon));
211 try expect(math.approxEqAbs(f32, pow(f32, 37.45, 3.3), 155736.703125, epsilon));
212 try expect(math.approxEqAbs(f32, pow(f32, 89.123, 3.3), 2722489.5, epsilon));
213 try expect(math.approxEqAbs(f32, pow(f32, -1.0, 1e10), 1.0, epsilon));
214
215 try expect(math.approxEqAbs(f64, pow(f64, 0.0, 3.3), 0.0, epsilon));
216 try expect(math.approxEqAbs(f64, pow(f64, 0.8923, 3.3), 0.686572, epsilon));
217 try expect(math.approxEqAbs(f64, pow(f64, 0.2, 3.3), 0.004936, epsilon));
218 try expect(math.approxEqAbs(f64, pow(f64, 1.5, 3.3), 3.811546, epsilon));
219 try expect(math.approxEqAbs(f64, pow(f64, 37.45, 3.3), 155736.7160616, epsilon));
220 try expect(math.approxEqAbs(f64, pow(f64, 89.123, 3.3), 2722490.231436, epsilon));
221 try expect(math.approxEqAbs(f64, pow(f64, -1.0, 1e20), 1.0, epsilon));
222}
223
224test "special" {
225 // pow(x, +-0) = 1 for any x
226 try expect(pow(f32, 4, 0.0) == 1.0);
227 try expect(pow(f32, 7, -0.0) == 1.0);
228 try expect(pow(f32, math.nan(f32), -0.0) == 1.0);
229 // pow(1, y) = 1 for any y
230 try expect(pow(f32, 1.0, 4) == 1.0);
231 try expect(pow(f32, 1.0, 7) == 1.0);
232 try expect(pow(f32, 1.0, -math.inf(f32)) == 1.0);
233 try expect(pow(f32, 1.0, math.nan(f32)) == 1.0);
234 // pow(x, 1) = x for any x
235 try expect(pow(f32, 45, 1.0) == 45);
236 try expect(pow(f32, -45, 1.0) == -45);
237 try expect(math.isPositiveZero(pow(f32, 0.0, 1.0)));
238 try expect(math.isNegativeZero(pow(f32, -0.0, 1.0)));
239 try expect(math.isPositiveInf(pow(f32, math.inf(f32), 1.0)));
240 try expect(math.isNan(pow(f32, math.nan(f32), 1.0)));
241 // pow(nan, y) = nan for any y != 0
242 try expect(math.isNan(pow(f32, math.nan(f32), 5.0)));
243 // pow(x, nan) = nan for any x != 1
244 try expect(math.isNan(pow(f32, 5.0, math.nan(f32))));
245 // pow(+-0, y) = +-inf for y an odd integer < 0
246 try expect(math.isPositiveInf(pow(f32, 0.0, -1.0)));
247 try expect(math.isNegativeInf(pow(f32, -0.0, -5.0)));
248 // pow(+-0, -inf) = +inf
249 try expect(math.isPositiveInf(pow(f32, 0.0, -math.inf(f32))));
250 try expect(math.isPositiveInf(pow(f32, -0.0, -math.inf(f32))));
251 // pow(+-0, +inf) = +0
252 try expect(math.isPositiveZero(pow(f32, 0.0, math.inf(f32))));
253 try expect(math.isPositiveZero(pow(f32, -0.0, math.inf(f32))));
254 // pow(+-0, y) = +inf for finite y < 0 and not an odd integer
255 try expect(math.isPositiveInf(pow(f32, 0.0, -2.0)));
256 try expect(math.isPositiveInf(pow(f32, -0.0, -2.0)));
257 try expect(math.isPositiveInf(pow(f32, 0.0, -5.2)));
258 try expect(math.isPositiveInf(pow(f32, -0.0, -0.5)));
259 // pow(+-0, y) = +-0 for finite y an odd integer > 0
260 try expect(math.isPositiveZero(pow(f32, 0.0, 3.0)));
261 try expect(math.isNegativeZero(pow(f32, -0.0, 5.0)));
262 // pow(+-0, y) = +0 for finite y > 0 and not an odd integer
263 try expect(math.isPositiveZero(pow(f32, 0.0, 2.0)));
264 try expect(math.isPositiveZero(pow(f32, -0.0, 2.0)));
265 try expect(math.isPositiveZero(pow(f32, 0.0, 5.2)));
266 try expect(math.isPositiveZero(pow(f32, -0.0, 0.5)));
267 // pow(-1, +-inf) = 1
268 try expect(pow(f32, -1.0, math.inf(f32)) == 1.0);
269 try expect(pow(f32, -1.0, -math.inf(f32)) == 1.0);
270 // pow(x, +inf) = +inf for |x| > 1
271 try expect(math.isPositiveInf(pow(f32, 1.2, math.inf(f32))));
272 try expect(math.isPositiveInf(pow(f32, -1.2, math.inf(f32))));
273 // pow(x, -inf) = +0 for |x| > 1
274 try expect(math.isPositiveZero(pow(f32, 1.2, -math.inf(f32))));
275 try expect(math.isPositiveZero(pow(f32, -1.2, -math.inf(f32))));
276 // pow(x, +inf) = +0 for |x| < 1
277 try expect(math.isPositiveZero(pow(f32, 0.2, math.inf(f32))));
278 try expect(math.isPositiveZero(pow(f32, -0.2, math.inf(f32))));
279 // pow(x, -inf) = +inf for |x| < 1
280 try expect(math.isPositiveInf(pow(f32, 0.2, -math.inf(f32))));
281 try expect(math.isPositiveInf(pow(f32, -0.2, -math.inf(f32))));
282 // pow(+inf, y) = +inf for y > 0
283 try expect(math.isPositiveInf(pow(f32, math.inf(f32), 2.0)));
284 try expect(math.isPositiveInf(pow(f32, math.inf(f32), 0.2)));
285 // pow(+inf, y) = +0 for y < 0
286 try expect(math.isPositiveZero(pow(f32, math.inf(f32), -2.0)));
287 try expect(math.isPositiveZero(pow(f32, math.inf(f32), -0.2)));
288 // pow(-inf, y) = -inf for y an odd integer > 0
289 try expect(math.isNegativeInf(pow(f32, -math.inf(f32), 5.0)));
290 // pow(-inf, +inf) = +inf
291 try expect(math.isPositiveInf(pow(f32, -math.inf(f32), math.inf(f32))));
292 // pow(-inf, -inf) = +0
293 try expect(math.isPositiveZero(pow(f32, -math.inf(f32), -math.inf(f32))));
294 // pow(-inf, y) = +inf for finite y > 0 and not an odd integer
295 try expect(math.isPositiveInf(pow(f32, -math.inf(f32), 4.0)));
296 try expect(math.isPositiveInf(pow(f32, -math.inf(f32), 0.5)));
297 // pow(-inf, y) = -0 for finite y an odd integer < 0
298 try expect(math.isNegativeZero(pow(f32, -math.inf(f32), -3.0)));
299 // pow(-inf, y) = +0 for finite y < 0 and not an odd integer
300 try expect(math.isPositiveZero(pow(f32, -math.inf(f32), -2.0)));
301 try expect(math.isPositiveZero(pow(f32, -math.inf(f32), -5.2)));
302 // pow(x, y) = nan for finite x < 0 and finite non-integer y
303 try expect(math.isNan(pow(f32, -1.0, 1.2)));
304 try expect(math.isNan(pow(f32, -12.4, -78.5)));
305}
306
307test "overflow" {
308 try expect(math.isPositiveInf(pow(f64, 2, 1 << 32)));
309 try expect(pow(f64, 2, -(1 << 32)) == 0);
310 try expect(math.isNegativeInf(pow(f64, -2, (1 << 32) + 1)));
311 try expect(pow(f64, 0.5, 1 << 45) == 0);
312 try expect(math.isPositiveInf(pow(f64, 0.5, -(1 << 45))));
313 try expect(math.isPositiveInf(pow(f64, -2, 1 << 64)));
314 try expect(pow(f64, 0.5, 1 << 64) == 0);
315}