authorgravatar for marc@tiehu.isMarc Tiehuis <marc@tiehu.is> 2021-11-23 20:10:32+13:00
committergravatar for andrew@ziglang.orgAndrew Kelley <andrew@ziglang.org> 2021-11-23 14:47:01-05:00
log8f1e417757f2172485cdea792580e2a2afb53ffa
tree42b7f53a6b872f5c8afaec97f24785bbb0e5d5cb
parent8e6c038dd81fa7eddfdac882aa3b1d5edbdcf329

std/math: add ldexp and make scalbn an alias

We assume we are compiled on a base-2 radix floating point system. This is a reasonable assumption. musl libc as an example also assumes this. We implement scalbn as an alias for ldexp, since ldexp is defined as 2 regardless of the float radix. This is opposite to musl which defines scalbn in terms of ldexp. Closes #9799.

3 files changed, 98 insertions(+), 82 deletions(-)

lib/std/math.zig+1
......@@ -241,6 +241,7 @@ pub const isNegativeInf = @import("math/isinf.zig").isNegativeInf;
241241pub const isNormal = @import("math/isnormal.zig").isNormal;
242242pub const signbit = @import("math/signbit.zig").signbit;
243243pub const scalbn = @import("math/scalbn.zig").scalbn;
244pub const ldexp = @import("math/ldexp.zig").ldexp;
244245pub const pow = @import("math/pow.zig").pow;
245246pub const powi = @import("math/powi.zig").powi;
246247pub const sqrt = @import("math/sqrt.zig").sqrt;
lib/std/math/ldexp.zig created+92
......@@ -0,0 +1,92 @@
1// Ported from musl, which is licensed under the MIT license:
2// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
3//
4// https://git.musl-libc.org/cgit/musl/tree/src/math/ldexpf.c
5// https://git.musl-libc.org/cgit/musl/tree/src/math/ldexp.c
6
7const std = @import("std");
8const math = std.math;
9const assert = std.debug.assert;
10const expect = std.testing.expect;
11
12/// Returns x * 2^n.
13pub fn ldexp(x: anytype, n: i32) @TypeOf(x) {
14 var base = x;
15 var shift = n;
16
17 const T = @TypeOf(base);
18 const IntT = std.meta.Int(.unsigned, @bitSizeOf(T));
19 if (@typeInfo(T) != .Float) {
20 @compileError("ldexp not implemented for " ++ @typeName(T));
21 }
22
23 const mantissa_bits = math.floatMantissaBits(T);
24 const exponent_bits = math.floatExponentBits(T);
25 const exponent_bias = (1 << (exponent_bits - 1)) - 1;
26 const exponent_min = 1 - exponent_bias;
27 const exponent_max = exponent_bias;
28
29 // fix double rounding errors in subnormal ranges
30 // https://git.musl-libc.org/cgit/musl/commit/src/math/ldexp.c?id=8c44a060243f04283ca68dad199aab90336141db
31 const scale_min_expo = exponent_min + mantissa_bits + 1;
32 const scale_min = @bitCast(T, @as(IntT, scale_min_expo + exponent_bias) << mantissa_bits);
33 const scale_max = @bitCast(T, @intCast(IntT, exponent_max + exponent_bias) << mantissa_bits);
34
35 // scale `shift` within floating point limits, if possible
36 // second pass is possible due to subnormal range
37 // third pass always results in +/-0.0 or +/-inf
38 if (shift > exponent_max) {
39 base *= scale_max;
40 shift -= exponent_max;
41 if (shift > exponent_max) {
42 base *= scale_max;
43 shift -= exponent_max;
44 if (shift > exponent_max) shift = exponent_max;
45 }
46 } else if (shift < exponent_min) {
47 base *= scale_min;
48 shift -= scale_min_expo;
49 if (shift < exponent_min) {
50 base *= scale_min;
51 shift -= scale_min_expo;
52 if (shift < exponent_min) shift = exponent_min;
53 }
54 }
55
56 return base * @bitCast(T, @intCast(IntT, shift + exponent_bias) << mantissa_bits);
57}
58
59test "math.ldexp" {
60 // basic usage
61 try expect(ldexp(@as(f16, 1.5), 4) == 24.0);
62 try expect(ldexp(@as(f32, 1.5), 4) == 24.0);
63 try expect(ldexp(@as(f64, 1.5), 4) == 24.0);
64 try expect(ldexp(@as(f128, 1.5), 4) == 24.0);
65
66 // subnormals
67 try expect(math.isNormal(ldexp(@as(f16, 1.0), -14)));
68 try expect(!math.isNormal(ldexp(@as(f16, 1.0), -15)));
69 try expect(math.isNormal(ldexp(@as(f32, 1.0), -126)));
70 try expect(!math.isNormal(ldexp(@as(f32, 1.0), -127)));
71 try expect(math.isNormal(ldexp(@as(f64, 1.0), -1022)));
72 try expect(!math.isNormal(ldexp(@as(f64, 1.0), -1023)));
73 try expect(math.isNormal(ldexp(@as(f128, 1.0), -16382)));
74 try expect(!math.isNormal(ldexp(@as(f128, 1.0), -16383)));
75 // unreliable due to lack of native f16 support, see talk on PR #8733
76 // try expect(ldexp(@as(f16, 0x1.1FFp-1), -14 - 9) == math.f16_true_min);
77 try expect(ldexp(@as(f32, 0x1.3FFFFFp-1), -126 - 22) == math.f32_true_min);
78 try expect(ldexp(@as(f64, 0x1.7FFFFFFFFFFFFp-1), -1022 - 51) == math.f64_true_min);
79 try expect(ldexp(@as(f128, 0x1.7FFFFFFFFFFFFFFFFFFFFFFFFFFFp-1), -16382 - 111) == math.f128_true_min);
80
81 // float limits
82 try expect(ldexp(@as(f32, math.f32_max), -128 - 149) > 0.0);
83 try expect(ldexp(@as(f32, math.f32_max), -128 - 149 - 1) == 0.0);
84 try expect(!math.isPositiveInf(ldexp(@as(f16, math.f16_true_min), 15 + 24)));
85 try expect(math.isPositiveInf(ldexp(@as(f16, math.f16_true_min), 15 + 24 + 1)));
86 try expect(!math.isPositiveInf(ldexp(@as(f32, math.f32_true_min), 127 + 149)));
87 try expect(math.isPositiveInf(ldexp(@as(f32, math.f32_true_min), 127 + 149 + 1)));
88 try expect(!math.isPositiveInf(ldexp(@as(f64, math.f64_true_min), 1023 + 1074)));
89 try expect(math.isPositiveInf(ldexp(@as(f64, math.f64_true_min), 1023 + 1074 + 1)));
90 try expect(!math.isPositiveInf(ldexp(@as(f128, math.f128_true_min), 16383 + 16494)));
91 try expect(math.isPositiveInf(ldexp(@as(f128, math.f128_true_min), 16383 + 16494 + 1)));
92}
lib/std/math/scalbn.zig+5-82
......@@ -1,92 +1,15 @@
1// Ported from musl, which is licensed under the MIT license:
2// https://git.musl-libc.org/cgit/musl/tree/COPYRIGHT
3//
4// https://git.musl-libc.org/cgit/musl/tree/src/math/scalbnf.c
5// https://git.musl-libc.org/cgit/musl/tree/src/math/scalbn.c
6
71const std = @import("std");
8const math = std.math;
9const assert = std.debug.assert;
102const expect = std.testing.expect;
113
12/// Returns x * 2^n.
13pub fn scalbn(x: anytype, n: i32) @TypeOf(x) {
14 var base = x;
15 var shift = n;
16
17 const T = @TypeOf(base);
18 const IntT = std.meta.Int(.unsigned, @bitSizeOf(T));
19 if (@typeInfo(T) != .Float) {
20 @compileError("scalbn not implemented for " ++ @typeName(T));
21 }
22
23 const mantissa_bits = math.floatMantissaBits(T);
24 const exponent_bits = math.floatExponentBits(T);
25 const exponent_bias = (1 << (exponent_bits - 1)) - 1;
26 const exponent_min = 1 - exponent_bias;
27 const exponent_max = exponent_bias;
28
29 // fix double rounding errors in subnormal ranges
30 // https://git.musl-libc.org/cgit/musl/commit/src/math/scalbn.c?id=8c44a060243f04283ca68dad199aab90336141db
31 const scale_min_expo = exponent_min + mantissa_bits + 1;
32 const scale_min = @bitCast(T, @as(IntT, scale_min_expo + exponent_bias) << mantissa_bits);
33 const scale_max = @bitCast(T, @intCast(IntT, exponent_max + exponent_bias) << mantissa_bits);
34
35 // scale `shift` within floating point limits, if possible
36 // second pass is possible due to subnormal range
37 // third pass always results in +/-0.0 or +/-inf
38 if (shift > exponent_max) {
39 base *= scale_max;
40 shift -= exponent_max;
41 if (shift > exponent_max) {
42 base *= scale_max;
43 shift -= exponent_max;
44 if (shift > exponent_max) shift = exponent_max;
45 }
46 } else if (shift < exponent_min) {
47 base *= scale_min;
48 shift -= scale_min_expo;
49 if (shift < exponent_min) {
50 base *= scale_min;
51 shift -= scale_min_expo;
52 if (shift < exponent_min) shift = exponent_min;
53 }
54 }
55
56 return base * @bitCast(T, @intCast(IntT, shift + exponent_bias) << mantissa_bits);
57}
4/// Returns a * FLT_RADIX ^ exp.
5///
6/// Zig only supports binary radix IEEE-754 floats. Hence FLT_RADIX=2, and this is an alias for ldexp.
7pub const scalbn = @import("ldexp.zig").ldexp;
588
599test "math.scalbn" {
60 // basic usage
10 // Verify we are using radix 2.
6111 try expect(scalbn(@as(f16, 1.5), 4) == 24.0);
6212 try expect(scalbn(@as(f32, 1.5), 4) == 24.0);
6313 try expect(scalbn(@as(f64, 1.5), 4) == 24.0);
6414 try expect(scalbn(@as(f128, 1.5), 4) == 24.0);
65
66 // subnormals
67 try expect(math.isNormal(scalbn(@as(f16, 1.0), -14)));
68 try expect(!math.isNormal(scalbn(@as(f16, 1.0), -15)));
69 try expect(math.isNormal(scalbn(@as(f32, 1.0), -126)));
70 try expect(!math.isNormal(scalbn(@as(f32, 1.0), -127)));
71 try expect(math.isNormal(scalbn(@as(f64, 1.0), -1022)));
72 try expect(!math.isNormal(scalbn(@as(f64, 1.0), -1023)));
73 try expect(math.isNormal(scalbn(@as(f128, 1.0), -16382)));
74 try expect(!math.isNormal(scalbn(@as(f128, 1.0), -16383)));
75 // unreliable due to lack of native f16 support, see talk on PR #8733
76 // try expect(scalbn(@as(f16, 0x1.1FFp-1), -14 - 9) == math.f16_true_min);
77 try expect(scalbn(@as(f32, 0x1.3FFFFFp-1), -126 - 22) == math.f32_true_min);
78 try expect(scalbn(@as(f64, 0x1.7FFFFFFFFFFFFp-1), -1022 - 51) == math.f64_true_min);
79 try expect(scalbn(@as(f128, 0x1.7FFFFFFFFFFFFFFFFFFFFFFFFFFFp-1), -16382 - 111) == math.f128_true_min);
80
81 // float limits
82 try expect(scalbn(@as(f32, math.f32_max), -128 - 149) > 0.0);
83 try expect(scalbn(@as(f32, math.f32_max), -128 - 149 - 1) == 0.0);
84 try expect(!math.isPositiveInf(scalbn(@as(f16, math.f16_true_min), 15 + 24)));
85 try expect(math.isPositiveInf(scalbn(@as(f16, math.f16_true_min), 15 + 24 + 1)));
86 try expect(!math.isPositiveInf(scalbn(@as(f32, math.f32_true_min), 127 + 149)));
87 try expect(math.isPositiveInf(scalbn(@as(f32, math.f32_true_min), 127 + 149 + 1)));
88 try expect(!math.isPositiveInf(scalbn(@as(f64, math.f64_true_min), 1023 + 1074)));
89 try expect(math.isPositiveInf(scalbn(@as(f64, math.f64_true_min), 1023 + 1074 + 1)));
90 try expect(!math.isPositiveInf(scalbn(@as(f128, math.f128_true_min), 16383 + 16494)));
91 try expect(math.isPositiveInf(scalbn(@as(f128, math.f128_true_min), 16383 + 16494 + 1)));
9215}