authorgravatar for viridianmeow@gmail.comviri <viridianmeow@gmail.com> 2021-05-14 14:15:53-06:00
committergravatar for noreply@github.comGitHub <noreply@github.com> 2021-05-14 16:15:53-04:00
log612ad779bb2df8e7e15ab20dd3df3a2900ff82f0
tree8d37d819238fb9f062225189b2093395e435fba9
parent114c6612cb0709a317f941b49daa4d4849dc24e1
signaturebadge-question-mark Signed by PGP key 4AEE18F83AFDEB23

std: rework math.scalbn (#8733)

* std: fix overflow in math.scalbn32 * std: rewrite math.scalbn to be generic * std: support f128 in math.isNormal * std: enable f128 tests in math.scalbn

2 files changed, 75 insertions(+), 68 deletions(-)

lib/std/math/isnormal.zig+9-2
......@@ -14,16 +14,20 @@ pub fn isNormal(x: anytype) bool {
1414 switch (T) {
1515 f16 => {
1616 const bits = @bitCast(u16, x);
17 return (bits + 1024) & 0x7FFF >= 2048;
17 return (bits + (1 << 10)) & (maxInt(u16) >> 1) >= (1 << 11);
1818 },
1919 f32 => {
2020 const bits = @bitCast(u32, x);
21 return (bits + 0x00800000) & 0x7FFFFFFF >= 0x01000000;
21 return (bits + (1 << 23)) & (maxInt(u32) >> 1) >= (1 << 24);
2222 },
2323 f64 => {
2424 const bits = @bitCast(u64, x);
2525 return (bits + (1 << 52)) & (maxInt(u64) >> 1) >= (1 << 53);
2626 },
27 f128 => {
28 const bits = @bitCast(u128, x);
29 return (bits + (1 << 112)) & (maxInt(u128) >> 1) >= (1 << 113);
30 },
2731 else => {
2832 @compileError("isNormal not implemented for " ++ @typeName(T));
2933 },
......@@ -34,10 +38,13 @@ test "math.isNormal" {
3438 try expect(!isNormal(math.nan(f16)));
3539 try expect(!isNormal(math.nan(f32)));
3640 try expect(!isNormal(math.nan(f64)));
41 try expect(!isNormal(math.nan(f128)));
3742 try expect(!isNormal(@as(f16, 0)));
3843 try expect(!isNormal(@as(f32, 0)));
3944 try expect(!isNormal(@as(f64, 0)));
45 try expect(!isNormal(@as(f128, 0)));
4046 try expect(isNormal(@as(f16, 1.0)));
4147 try expect(isNormal(@as(f32, 1.0)));
4248 try expect(isNormal(@as(f64, 1.0)));
49 try expect(isNormal(@as(f128, 1.0)));
4350}
lib/std/math/scalbn.zig+66-66
......@@ -9,89 +9,89 @@
99// https://git.musl-libc.org/cgit/musl/tree/src/math/scalbnf.c
1010// https://git.musl-libc.org/cgit/musl/tree/src/math/scalbn.c
1111
12const std = @import("../std.zig");
12const std = @import("std");
1313const math = std.math;
14const assert = std.debug.assert;
1415const expect = std.testing.expect;
1516
1617/// Returns x * 2^n.
1718pub fn scalbn(x: anytype, n: i32) @TypeOf(x) {
18 const T = @TypeOf(x);
19 return switch (T) {
20 f32 => scalbn32(x, n),
21 f64 => scalbn64(x, n),
22 else => @compileError("scalbn not implemented for " ++ @typeName(T)),
23 };
24}
25
26fn scalbn32(x: f32, n_: i32) f32 {
27 var y = x;
28 var n = n_;
19 var base = x;
20 var shift = n;
2921
30 if (n > 127) {
31 y *= 0x1.0p127;
32 n -= 127;
33 if (n > 1023) {
34 y *= 0x1.0p127;
35 n -= 127;
36 if (n > 127) {
37 n = 127;
38 }
39 }
40 } else if (n < -126) {
41 y *= 0x1.0p-126 * 0x1.0p24;
42 n += 126 - 24;
43 if (n < -126) {
44 y *= 0x1.0p-126 * 0x1.0p24;
45 n += 126 - 24;
46 if (n < -126) {
47 n = -126;
48 }
49 }
22 const T = @TypeOf(base);
23 const IntT = std.meta.Int(.unsigned, @bitSizeOf(T));
24 if (@typeInfo(T) != .Float) {
25 @compileError("scalbn not implemented for " ++ @typeName(T));
5026 }
5127
52 const u = @intCast(u32, n +% 0x7F) << 23;
53 return y * @bitCast(f32, u);
54}
28 const mantissa_bits = math.floatMantissaBits(T);
29 const exponent_bits = math.floatExponentBits(T);
30 const exponent_bias = (1 << (exponent_bits - 1)) - 1;
31 const exponent_min = 1 - exponent_bias;
32 const exponent_max = exponent_bias;
5533
56fn scalbn64(x: f64, n_: i32) f64 {
57 var y = x;
58 var n = n_;
34 // fix double rounding errors in subnormal ranges
35 // https://git.musl-libc.org/cgit/musl/commit/src/math/scalbn.c?id=8c44a060243f04283ca68dad199aab90336141db
36 const scale_min_expo = exponent_min + mantissa_bits + 1;
37 const scale_min = @bitCast(T, @as(IntT, scale_min_expo + exponent_bias) << mantissa_bits);
38 const scale_max = @bitCast(T, @intCast(IntT, exponent_max + exponent_bias) << mantissa_bits);
5939
60 if (n > 1023) {
61 y *= 0x1.0p1023;
62 n -= 1023;
63 if (n > 1023) {
64 y *= 0x1.0p1023;
65 n -= 1023;
66 if (n > 1023) {
67 n = 1023;
68 }
40 // scale `shift` within floating point limits, if possible
41 // second pass is possible due to subnormal range
42 // third pass always results in +/-0.0 or +/-inf
43 if (shift > exponent_max) {
44 base *= scale_max;
45 shift -= exponent_max;
46 if (shift > exponent_max) {
47 base *= scale_max;
48 shift -= exponent_max;
49 if (shift > exponent_max) shift = exponent_max;
6950 }
70 } else if (n < -1022) {
71 y *= 0x1.0p-1022 * 0x1.0p53;
72 n += 1022 - 53;
73 if (n < -1022) {
74 y *= 0x1.0p-1022 * 0x1.0p53;
75 n += 1022 - 53;
76 if (n < -1022) {
77 n = -1022;
78 }
51 } else if (shift < exponent_min) {
52 base *= scale_min;
53 shift -= scale_min_expo;
54 if (shift < exponent_min) {
55 base *= scale_min;
56 shift -= scale_min_expo;
57 if (shift < exponent_min) shift = exponent_min;
7958 }
8059 }
8160
82 const u = @intCast(u64, n +% 0x3FF) << 52;
83 return y * @bitCast(f64, u);
61 return base * @bitCast(T, @intCast(IntT, shift + exponent_bias) << mantissa_bits);
8462}
8563
8664test "math.scalbn" {
87 try expect(scalbn(@as(f32, 1.5), 4) == scalbn32(1.5, 4));
88 try expect(scalbn(@as(f64, 1.5), 4) == scalbn64(1.5, 4));
89}
65 // basic usage
66 try expect(scalbn(@as(f16, 1.5), 4) == 24.0);
67 try expect(scalbn(@as(f32, 1.5), 4) == 24.0);
68 try expect(scalbn(@as(f64, 1.5), 4) == 24.0);
69 try expect(scalbn(@as(f128, 1.5), 4) == 24.0);
9070
91test "math.scalbn32" {
92 try expect(scalbn32(1.5, 4) == 24.0);
93}
71 // subnormals
72 try expect(math.isNormal(scalbn(@as(f16, 1.0), -14)));
73 try expect(!math.isNormal(scalbn(@as(f16, 1.0), -15)));
74 try expect(math.isNormal(scalbn(@as(f32, 1.0), -126)));
75 try expect(!math.isNormal(scalbn(@as(f32, 1.0), -127)));
76 try expect(math.isNormal(scalbn(@as(f64, 1.0), -1022)));
77 try expect(!math.isNormal(scalbn(@as(f64, 1.0), -1023)));
78 try expect(math.isNormal(scalbn(@as(f128, 1.0), -16382)));
79 try expect(!math.isNormal(scalbn(@as(f128, 1.0), -16383)));
80 // unreliable due to lack of native f16 support, see talk on PR #8733
81 // try expect(scalbn(@as(f16, 0x1.1FFp-1), -14 - 9) == math.f16_true_min);
82 try expect(scalbn(@as(f32, 0x1.3FFFFFp-1), -126 - 22) == math.f32_true_min);
83 try expect(scalbn(@as(f64, 0x1.7FFFFFFFFFFFFp-1), -1022 - 51) == math.f64_true_min);
84 try expect(scalbn(@as(f128, 0x1.7FFFFFFFFFFFFFFFFFFFFFFFFFFFp-1), -16382 - 111) == math.f128_true_min);
9485
95test "math.scalbn64" {
96 try expect(scalbn64(1.5, 4) == 24.0);
86 // float limits
87 try expect(scalbn(@as(f32, math.f32_max), -128 - 149) > 0.0);
88 try expect(scalbn(@as(f32, math.f32_max), -128 - 149 - 1) == 0.0);
89 try expect(!math.isPositiveInf(scalbn(@as(f16, math.f16_true_min), 15 + 24)));
90 try expect(math.isPositiveInf(scalbn(@as(f16, math.f16_true_min), 15 + 24 + 1)));
91 try expect(!math.isPositiveInf(scalbn(@as(f32, math.f32_true_min), 127 + 149)));
92 try expect(math.isPositiveInf(scalbn(@as(f32, math.f32_true_min), 127 + 149 + 1)));
93 try expect(!math.isPositiveInf(scalbn(@as(f64, math.f64_true_min), 1023 + 1074)));
94 try expect(math.isPositiveInf(scalbn(@as(f64, math.f64_true_min), 1023 + 1074 + 1)));
95 try expect(!math.isPositiveInf(scalbn(@as(f128, math.f128_true_min), 16383 + 16494)));
96 try expect(math.isPositiveInf(scalbn(@as(f128, math.f128_true_min), 16383 + 16494 + 1)));
9797}