Skip to content

Commit a2df84d

Browse files
SeanTUTSeanTUT
authored
std.math: rework modf
- Changed `modf_result` to `Modf` to better fit naming conventions - Reworked `modf` to be far simpler and support all floating point types (as well as vectors) (I have done benchmarks and can confirm that the performance is roughly equivalent to the old implementation) - Added more descriptive tests for modf - Deprecated `modf32_result` and `modf64_result` in favor of `Modf(f32)` and `Modf(f64)` respectively
1 parent 2d443cd commit a2df84d

File tree

2 files changed

+103
-169
lines changed

2 files changed

+103
-169
lines changed

‎lib/std/math.zig‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -112,6 +112,8 @@ pub const qnan_f80 = @compileError("Deprecated: use `nan(f80)` instead");
112112
pub const qnan_u128 = @compileError("Deprecated: use `@as(u128, @bitCast(nan(f128)))` instead");
113113
pub const qnan_f128 = @compileError("Deprecated: use `nan(f128)` instead");
114114
pub const epsilon = @compileError("Deprecated: use `floatEps` instead");
115+
pub const modf32_result = @compileError("Deprecated: use `Modf(f32)` instead");
116+
pub const modf64_result = @compileError("Deprecated: use `Modf(f64)` instead");
115117

116118
/// Performs an approximate comparison of two floating point values `x` and `y`.
117119
/// Returns true if the absolute difference between them is less or equal than
@@ -255,8 +257,7 @@ pub const isSignalNan = @import("math/isnan.zig").isSignalNan;
255257
pub const frexp = @import("math/frexp.zig").frexp;
256258
pub const Frexp = @import("math/frexp.zig").Frexp;
257259
pub const modf = @import("math/modf.zig").modf;
258-
pub const modf32_result = @import("math/modf.zig").modf32_result;
259-
pub const modf64_result = @import("math/modf.zig").modf64_result;
260+
pub const Modf = @import("math/modf.zig").Modf;
260261
pub const copysign = @import("math/copysign.zig").copysign;
261262
pub const isFinite = @import("math/isfinite.zig").isFinite;
262263
pub const isInf = @import("math/isinf.zig").isInf;
@@ -418,8 +419,7 @@ test {
418419
_ = frexp;
419420
_ = Frexp;
420421
_ = modf;
421-
_ = modf32_result;
422-
_ = modf64_result;
422+
_ = Modf;
423423
_ = copysign;
424424
_ = isFinite;
425425
_ = isInf;

‎lib/std/math/modf.zig‎

Lines changed: 99 additions & 165 deletions
Original file line numberDiff line numberDiff line change
@@ -1,207 +1,141 @@
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/modff.c
5-
// https://git.musl-libc.org/cgit/musl/tree/src/math/modf.c
6-
71
const std = @import("../std.zig");
82
const math = std.math;
93
const expect = std.testing.expect;
104
const expectEqual = std.testing.expectEqual;
11-
const maxInt = std.math.maxInt;
5+
const expectApproxEqAbs = std.testing.expectApproxEqAbs;
126

13-
fn modf_result(comptime T: type) type {
7+
pub fn Modf(comptime T: type) type {
148
return struct {
159
fpart: T,
1610
ipart: T,
1711
};
1812
}
19-
pub const modf32_result = modf_result(f32);
20-
pub const modf64_result = modf_result(f64);
2113

2214
/// Returns the integer and fractional floating-point numbers that sum to x. The sign of each
2315
/// result is the same as the sign of x.
16+
/// In comptime, may be used with comptime_float
2417
///
2518
/// Special Cases:
2619
/// - modf(+-inf) = +-inf, nan
2720
/// - modf(nan) = nan, nan
28-
pub fn modf(x: anytype) modf_result(@TypeOf(x)) {
29-
const T = @TypeOf(x);
30-
return switch (T) {
31-
f32 => modf32(x),
32-
f64 => modf64(x),
33-
else => @compileError("modf not implemented for " ++ @typeName(T)),
21+
pub fn modf(x: anytype) Modf(@TypeOf(x)) {
22+
const ipart = @trunc(x);
23+
return .{
24+
.ipart = ipart,
25+
.fpart = x - ipart,
3426
};
3527
}
3628

37-
fn modf32(x: f32) modf32_result {
38-
var result: modf32_result = undefined;
29+
test modf {
30+
inline for ([_]type{ f16, f32, f64, f80, f128 }) |T| {
31+
const epsilon: comptime_float = @max(1e-6, math.floatEps(T));
3932

40-
const u: u32 = @bitCast(x);
41-
const e = @as(i32, @intCast((u >> 23) & 0xFF)) - 0x7F;
42-
const us = u & 0x80000000;
33+
var r: Modf(T) = undefined;
4334

44-
// TODO: Shouldn't need this.
45-
if (math.isInf(x)) {
46-
result.ipart = x;
47-
result.fpart = math.nan(f32);
48-
return result;
49-
}
35+
r = modf(@as(T, 1.0));
36+
try expectEqual(1.0, r.ipart);
37+
try expectEqual(0.0, r.fpart);
5038

51-
// no fractional part
52-
if (e >= 23) {
53-
result.ipart = x;
54-
if (e == 0x80 and u << 9 != 0) { // nan
55-
result.fpart = x;
56-
} else {
57-
result.fpart = @as(f32, @bitCast(us));
58-
}
59-
return result;
60-
}
39+
r = modf(@as(T, 0.34682));
40+
try expectEqual(0.0, r.ipart);
41+
try expectApproxEqAbs(@as(T, 0.34682), r.fpart, epsilon);
6142

62-
// no integral part
63-
if (e < 0) {
64-
result.ipart = @as(f32, @bitCast(us));
65-
result.fpart = x;
66-
return result;
67-
}
43+
r = modf(@as(T, 2.54576));
44+
try expectEqual(2.0, r.ipart);
45+
try expectApproxEqAbs(0.54576, r.fpart, epsilon);
6846

69-
const mask = @as(u32, 0x007FFFFF) >> @as(u5, @intCast(e));
70-
if (u & mask == 0) {
71-
result.ipart = x;
72-
result.fpart = @as(f32, @bitCast(us));
73-
return result;
47+
r = modf(@as(T, 3.9782));
48+
try expectEqual(3.0, r.ipart);
49+
try expectApproxEqAbs(0.9782, r.fpart, epsilon);
7450
}
75-
76-
const uf: f32 = @bitCast(u & ~mask);
77-
result.ipart = uf;
78-
result.fpart = x - uf;
79-
return result;
8051
}
8152

82-
fn modf64(x: f64) modf64_result {
83-
var result: modf64_result = undefined;
84-
85-
const u: u64 = @bitCast(x);
86-
const e = @as(i32, @intCast((u >> 52) & 0x7FF)) - 0x3FF;
87-
const us = u & (1 << 63);
88-
89-
if (math.isInf(x)) {
90-
result.ipart = x;
91-
result.fpart = math.nan(f64);
92-
return result;
93-
}
94-
95-
// no fractional part
96-
if (e >= 52) {
97-
result.ipart = x;
98-
if (e == 0x400 and u << 12 != 0) { // nan
99-
result.fpart = x;
100-
} else {
101-
result.fpart = @as(f64, @bitCast(us));
53+
/// Generate a namespace of tests for modf on values of the given type
54+
fn ModfTests(comptime T: type) type {
55+
return struct {
56+
test "normal" {
57+
const epsilon: comptime_float = @max(1e-6, math.floatEps(T));
58+
var r: Modf(T) = undefined;
59+
60+
r = modf(@as(T, 1.0));
61+
try expectEqual(1.0, r.ipart);
62+
try expectEqual(0.0, r.fpart);
63+
64+
r = modf(@as(T, 0.34682));
65+
try expectEqual(0.0, r.ipart);
66+
try expectApproxEqAbs(0.34682, r.fpart, epsilon);
67+
68+
r = modf(@as(T, 3.97812));
69+
try expectEqual(3.0, r.ipart);
70+
// account for precision error
71+
const expected_a: T = 3.97812 - @as(T, 3);
72+
try expectApproxEqAbs(expected_a, r.fpart, epsilon);
73+
74+
r = modf(@as(T, 43874.3));
75+
try expectEqual(43874.0, r.ipart);
76+
// account for precision error
77+
const expected_b: T = 43874.3 - @as(T, 43874);
78+
try expectApproxEqAbs(expected_b, r.fpart, epsilon);
79+
80+
r = modf(@as(T, 1234.340780));
81+
try expectEqual(1234.0, r.ipart);
82+
// account for precision error
83+
const expected_c: T = 1234.340780 - @as(T, 1234);
84+
try expectApproxEqAbs(expected_c, r.fpart, epsilon);
10285
}
103-
return result;
104-
}
105-
106-
// no integral part
107-
if (e < 0) {
108-
result.ipart = @as(f64, @bitCast(us));
109-
result.fpart = x;
110-
return result;
111-
}
112-
113-
const mask = @as(u64, maxInt(u64) >> 12) >> @as(u6, @intCast(e));
114-
if (u & mask == 0) {
115-
result.ipart = x;
116-
result.fpart = @as(f64, @bitCast(us));
117-
return result;
118-
}
119-
120-
const uf = @as(f64, @bitCast(u & ~mask));
121-
result.ipart = uf;
122-
result.fpart = x - uf;
123-
return result;
124-
}
86+
test "vector" {
87+
// Currently, a compiler bug is breaking the usage
88+
// of @trunc on @Vector types
12589

126-
test modf {
127-
const a = modf(@as(f32, 1.0));
128-
const b = modf32(1.0);
129-
// NOTE: No struct comparison on generic return type function? non-named, makes sense, but still.
130-
try expectEqual(a, b);
131-
}
132-
133-
test modf32 {
134-
const epsilon = 0.000001;
135-
var r: modf32_result = undefined;
136-
137-
r = modf32(1.0);
138-
try expect(math.approxEqAbs(f32, r.ipart, 1.0, epsilon));
139-
try expect(math.approxEqAbs(f32, r.fpart, 0.0, epsilon));
140-
141-
r = modf32(2.545);
142-
try expect(math.approxEqAbs(f32, r.ipart, 2.0, epsilon));
143-
try expect(math.approxEqAbs(f32, r.fpart, 0.545, epsilon));
90+
// TODO: Repopulate the below array and
91+
// remove the skip statement once this
92+
// bug is fixed
14493

145-
r = modf32(3.978123);
146-
try expect(math.approxEqAbs(f32, r.ipart, 3.0, epsilon));
147-
try expect(math.approxEqAbs(f32, r.fpart, 0.978123, epsilon));
94+
// const widths = [_]comptime_int{ 1, 2, 3, 4, 8, 16 };
95+
const widths = [_]comptime_int{};
14896

149-
r = modf32(43874.3);
150-
try expect(math.approxEqAbs(f32, r.ipart, 43874, epsilon));
151-
try expect(math.approxEqAbs(f32, r.fpart, 0.300781, epsilon));
97+
if (widths.len == 0)
98+
return error.SkipZigTest;
15299

153-
r = modf32(1234.340780);
154-
try expect(math.approxEqAbs(f32, r.ipart, 1234, epsilon));
155-
try expect(math.approxEqAbs(f32, r.fpart, 0.340820, epsilon));
156-
}
157-
158-
test modf64 {
159-
const epsilon = 0.000001;
160-
var r: modf64_result = undefined;
161-
162-
r = modf64(1.0);
163-
try expect(math.approxEqAbs(f64, r.ipart, 1.0, epsilon));
164-
try expect(math.approxEqAbs(f64, r.fpart, 0.0, epsilon));
100+
inline for (widths) |len| {
101+
const V: type = @Vector(len, T);
102+
var r: Modf(V) = undefined;
165103

166-
r = modf64(2.545);
167-
try expect(math.approxEqAbs(f64, r.ipart, 2.0, epsilon));
168-
try expect(math.approxEqAbs(f64, r.fpart, 0.545, epsilon));
104+
r = modf(@as(V, @splat(1.0)));
105+
try expectEqual(@as(V, @splat(1.0)), r.ipart);
106+
try expectEqual(@as(V, @splat(0.0)), r.fpart);
169107

170-
r = modf64(3.978123);
171-
try expect(math.approxEqAbs(f64, r.ipart, 3.0, epsilon));
172-
try expect(math.approxEqAbs(f64, r.fpart, 0.978123, epsilon));
108+
r = modf(@as(V, @splat(2.75)));
109+
try expectEqual(@as(V, @splat(2.0)), r.ipart);
110+
try expectEqual(@as(V, @splat(0.75)), r.fpart);
173111

174-
r = modf64(43874.3);
175-
try expect(math.approxEqAbs(f64, r.ipart, 43874, epsilon));
176-
try expect(math.approxEqAbs(f64, r.fpart, 0.3, epsilon));
177-
178-
r = modf64(1234.340780);
179-
try expect(math.approxEqAbs(f64, r.ipart, 1234, epsilon));
180-
try expect(math.approxEqAbs(f64, r.fpart, 0.340780, epsilon));
181-
}
112+
r = modf(@as(V, @splat(0.2)));
113+
try expectEqual(@as(V, @splat(0.0)), r.ipart);
114+
try expectEqual(@as(V, @splat(0.2)), r.fpart);
182115

183-
test "modf32.special" {
184-
var r: modf32_result = undefined;
185-
186-
r = modf32(math.inf(f32));
187-
try expect(math.isPositiveInf(r.ipart) and math.isNan(r.fpart));
116+
r = modf(std.simd.iota(T, len) + @as(V, @splat(0.5)));
117+
try expectEqual(std.simd.iota(T, len), r.ipart);
118+
try expectEqual(@as(V, @splat(0.5)), r.fpart);
119+
}
120+
}
121+
test "inf" {
122+
var r: Modf(T) = undefined;
188123

189-
r = modf32(-math.inf(f32));
190-
try expect(math.isNegativeInf(r.ipart) and math.isNan(r.fpart));
124+
r = modf(math.inf(T));
125+
try expect(math.isPositiveInf(r.ipart) and math.isNan(r.fpart));
191126

192-
r = modf32(math.nan(f32));
193-
try expect(math.isNan(r.ipart) and math.isNan(r.fpart));
127+
r = modf(-math.inf(T));
128+
try expect(math.isNegativeInf(r.ipart) and math.isNan(r.fpart));
129+
}
130+
test "nan" {
131+
const r: Modf(T) = modf(math.nan(T));
132+
try expect(math.isNan(r.ipart) and math.isNan(r.fpart));
133+
}
134+
};
194135
}
195136

196-
test "modf64.special" {
197-
var r: modf64_result = undefined;
198-
199-
r = modf64(math.inf(f64));
200-
try expect(math.isPositiveInf(r.ipart) and math.isNan(r.fpart));
201-
202-
r = modf64(-math.inf(f64));
203-
try expect(math.isNegativeInf(r.ipart) and math.isNan(r.fpart));
204-
205-
r = modf64(math.nan(f64));
206-
try expect(math.isNan(r.ipart) and math.isNan(r.fpart));
137+
comptime {
138+
for ([_]type{ f16, f32, f64, f80, f128 }) |T| {
139+
_ = ModfTests(T);
140+
}
207141
}

0 commit comments

Comments
 (0)