Skip to content

Commit 8f1e417

Browse files
tiehuisandrewrk
tiehuis
authored and
andrewrk
committed
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.
1 parent 8e6c038 commit 8f1e417

File tree

3 files changed

+98
-82
lines changed

3 files changed

+98
-82
lines changed

‎lib/std/math.zig

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -241,6 +241,7 @@ pub const isNegativeInf = @import("math/isinf.zig").isNegativeInf;
241241
pub const isNormal = @import("math/isnormal.zig").isNormal;
242242
pub const signbit = @import("math/signbit.zig").signbit;
243243
pub const scalbn = @import("math/scalbn.zig").scalbn;
244+
pub const ldexp = @import("math/ldexp.zig").ldexp;
244245
pub const pow = @import("math/pow.zig").pow;
245246
pub const powi = @import("math/powi.zig").powi;
246247
pub const sqrt = @import("math/sqrt.zig").sqrt;

‎lib/std/math/ldexp.zig

Lines changed: 92 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -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+
7+
const std = @import("std");
8+
const math = std.math;
9+
const assert = std.debug.assert;
10+
const expect = std.testing.expect;
11+
12+
/// Returns x * 2^n.
13+
pub 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+
59+
test "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

Lines changed: 5 additions & 82 deletions
Original file line numberDiff line numberDiff line change
@@ -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-
71
const std = @import("std");
8-
const math = std.math;
9-
const assert = std.debug.assert;
102
const expect = std.testing.expect;
113

12-
/// Returns x * 2^n.
13-
pub 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.
7+
pub const scalbn = @import("ldexp.zig").ldexp;
588

599
test "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
}

0 commit comments

Comments
 (0)