| ... | ... | @@ -13,13 +13,13 @@ const expect = std.testing.expect; |
| 13 | 13 | /// - pow(x, +-0) = 1 for any x |
| 14 | 14 | /// - pow(1, y) = 1 for any y |
| 15 | 15 | /// - pow(x, 1) = x for any x |
| 16 | | /// - pow(nan, y) = nan |
| 17 | | /// - pow(x, nan) = nan |
| 16 | /// - pow(nan, y) = nan for any y != 0 |
| 17 | /// - pow(x, nan) = nan for any x != 1 |
| 18 | 18 | /// - pow(+-0, y) = +-inf for y an odd integer < 0 |
| 19 | 19 | /// - pow(+-0, -inf) = +inf |
| 20 | 20 | /// - pow(+-0, +inf) = +0 |
| 21 | 21 | /// - pow(+-0, y) = +inf for finite y < 0 and not an odd integer |
| 22 | | /// - pow(+-0, y) = +-0 for y an odd integer > 0 |
| 22 | /// - pow(+-0, y) = +-0 for finite y an odd integer > 0 |
| 23 | 23 | /// - pow(+-0, y) = +0 for finite y > 0 and not an odd integer |
| 24 | 24 | /// - pow(-1, +-inf) = 1 |
| 25 | 25 | /// - pow(x, +inf) = +inf for |x| > 1 |
| ... | ... | @@ -39,20 +39,20 @@ pub fn pow(comptime T: type, x: T, y: T) T { |
| 39 | 39 | @compileError("pow not implemented for " ++ @typeName(T)); |
| 40 | 40 | } |
| 41 | 41 | |
| 42 | | // pow(x, +-0) = 1 for all x |
| 43 | | // pow(1, y) = 1 for all y |
| 42 | // pow(x, +-0) = 1 for any x |
| 43 | // pow(1, y) = 1 for any y |
| 44 | 44 | if (y == 0 or x == 1) { |
| 45 | 45 | return 1; |
| 46 | 46 | } |
| 47 | 47 | |
| 48 | | // pow(nan, y) = nan for all y |
| 49 | | // pow(x, nan) = nan for all x |
| 48 | // pow(nan, y) = nan for any y != 0 |
| 49 | // pow(x, nan) = nan for any x != 1 |
| 50 | 50 | if (math.isNan(x) or math.isNan(y)) { |
| 51 | 51 | @branchHint(.unlikely); |
| 52 | 52 | return math.nan(T); |
| 53 | 53 | } |
| 54 | 54 | |
| 55 | | // pow(x, 1) = x for all x |
| 55 | // pow(x, 1) = x for any x |
| 56 | 56 | if (y == 1) { |
| 57 | 57 | return x; |
| 58 | 58 | } |
| ... | ... | @@ -77,7 +77,7 @@ pub fn pow(comptime T: type, x: T, y: T) T { |
| 77 | 77 | } |
| 78 | 78 | |
| 79 | 79 | if (math.isInf(y)) { |
| 80 | | // pow(-1, inf) = 1 for all x |
| 80 | // pow(-1, inf) = 1 |
| 81 | 81 | if (x == -1) { |
| 82 | 82 | return 1.0; |
| 83 | 83 | } |
| ... | ... | @@ -124,7 +124,9 @@ pub fn pow(comptime T: type, x: T, y: T) T { |
| 124 | 124 | return math.nan(T); |
| 125 | 125 | } |
| 126 | 126 | if (yi >= 1 << (@typeInfo(T).float.bits - 1)) { |
| 127 | | return @exp(y * @log(x)); |
| 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))); |
| 128 | 130 | } |
| 129 | 131 | |
| 130 | 132 | // a = a1 * 2^ae |
| ... | ... | @@ -208,6 +210,7 @@ test pow { |
| 208 | 210 | try expect(math.approxEqAbs(f32, pow(f32, 1.5, 3.3), 3.811546, epsilon)); |
| 209 | 211 | try expect(math.approxEqAbs(f32, pow(f32, 37.45, 3.3), 155736.703125, epsilon)); |
| 210 | 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)); |
| 211 | 214 | |
| 212 | 215 | try expect(math.approxEqAbs(f64, pow(f64, 0.0, 3.3), 0.0, epsilon)); |
| 213 | 216 | try expect(math.approxEqAbs(f64, pow(f64, 0.8923, 3.3), 0.686572, epsilon)); |
| ... | ... | @@ -215,48 +218,90 @@ test pow { |
| 215 | 218 | try expect(math.approxEqAbs(f64, pow(f64, 1.5, 3.3), 3.811546, epsilon)); |
| 216 | 219 | try expect(math.approxEqAbs(f64, pow(f64, 37.45, 3.3), 155736.7160616, epsilon)); |
| 217 | 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)); |
| 218 | 222 | } |
| 219 | 223 | |
| 220 | 224 | test "special" { |
| 221 | | const epsilon = 0.000001; |
| 222 | | |
| 225 | // pow(x, +-0) = 1 for any x |
| 223 | 226 | try expect(pow(f32, 4, 0.0) == 1.0); |
| 224 | 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 |
| 225 | 235 | try expect(pow(f32, 45, 1.0) == 45); |
| 226 | 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 |
| 227 | 242 | try expect(math.isNan(pow(f32, math.nan(f32), 5.0))); |
| 228 | | try expect(math.isPositiveInf(pow(f32, -math.inf(f32), 0.5))); |
| 229 | | try expect(math.isPositiveInf(pow(f32, -0.0, -0.5))); |
| 230 | | try expect(pow(f32, -0.0, 0.5) == 0); |
| 243 | // pow(x, nan) = nan for any x != 1 |
| 231 | 244 | try expect(math.isNan(pow(f32, 5.0, math.nan(f32)))); |
| 245 | // pow(+-0, y) = +-inf for y an odd integer < 0 |
| 232 | 246 | try expect(math.isPositiveInf(pow(f32, 0.0, -1.0))); |
| 233 | | //expect(math.isNegativeInf(pow(f32, -0.0, -3.0))); TODO is this required? |
| 247 | try expect(math.isNegativeInf(pow(f32, -0.0, -5.0))); |
| 248 | // pow(+-0, -inf) = +inf |
| 234 | 249 | try expect(math.isPositiveInf(pow(f32, 0.0, -math.inf(f32)))); |
| 235 | 250 | try expect(math.isPositiveInf(pow(f32, -0.0, -math.inf(f32)))); |
| 236 | | try expect(pow(f32, 0.0, math.inf(f32)) == 0.0); |
| 237 | | try expect(pow(f32, -0.0, math.inf(f32)) == 0.0); |
| 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 |
| 238 | 255 | try expect(math.isPositiveInf(pow(f32, 0.0, -2.0))); |
| 239 | 256 | try expect(math.isPositiveInf(pow(f32, -0.0, -2.0))); |
| 240 | | try expect(pow(f32, 0.0, 1.0) == 0.0); |
| 241 | | try expect(pow(f32, -0.0, 1.0) == -0.0); |
| 242 | | try expect(pow(f32, 0.0, 2.0) == 0.0); |
| 243 | | try expect(pow(f32, -0.0, 2.0) == 0.0); |
| 244 | | try expect(math.approxEqAbs(f32, pow(f32, -1.0, math.inf(f32)), 1.0, epsilon)); |
| 245 | | try expect(math.approxEqAbs(f32, pow(f32, -1.0, -math.inf(f32)), 1.0, epsilon)); |
| 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 |
| 246 | 271 | try expect(math.isPositiveInf(pow(f32, 1.2, math.inf(f32)))); |
| 247 | 272 | try expect(math.isPositiveInf(pow(f32, -1.2, math.inf(f32)))); |
| 248 | | try expect(pow(f32, 1.2, -math.inf(f32)) == 0.0); |
| 249 | | try expect(pow(f32, -1.2, -math.inf(f32)) == 0.0); |
| 250 | | try expect(pow(f32, 0.2, math.inf(f32)) == 0.0); |
| 251 | | try expect(pow(f32, -0.2, math.inf(f32)) == 0.0); |
| 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 |
| 252 | 280 | try expect(math.isPositiveInf(pow(f32, 0.2, -math.inf(f32)))); |
| 253 | 281 | try expect(math.isPositiveInf(pow(f32, -0.2, -math.inf(f32)))); |
| 254 | | try expect(math.isPositiveInf(pow(f32, math.inf(f32), 1.0))); |
| 255 | | try expect(pow(f32, math.inf(f32), -1.0) == 0.0); |
| 256 | | //expect(pow(f32, -math.inf(f32), 5.0) == pow(f32, -0.0, -5.0)); TODO support negative 0? |
| 257 | | try expect(pow(f32, -math.inf(f32), -5.2) == pow(f32, -0.0, 5.2)); |
| 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 |
| 258 | 303 | try expect(math.isNan(pow(f32, -1.0, 1.2))); |
| 259 | | try expect(math.isNan(pow(f32, -12.4, 78.5))); |
| 304 | try expect(math.isNan(pow(f32, -12.4, -78.5))); |
| 260 | 305 | } |
| 261 | 306 | |
| 262 | 307 | test "overflow" { |
| ... | ... | @@ -265,4 +310,6 @@ test "overflow" { |
| 265 | 310 | try expect(math.isNegativeInf(pow(f64, -2, (1 << 32) + 1))); |
| 266 | 311 | try expect(pow(f64, 0.5, 1 << 45) == 0); |
| 267 | 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); |
| 268 | 315 | } |