1const std = @import("../std.zig");
2const builtin = @import("builtin");
3const assert = std.debug.assert;
4const expect = std.testing.expect;
5const expectEqual = std.testing.expectEqual;
6
7/// A namespace for functions that deal with floats that provide greater than
8/// double precision (`f80`, `f128`, `c_longdouble`). Commonly referred to as
9/// `long double` in C.
10pub const long_double = struct {
11 const U80 = @Int(.unsigned, 80);
12
13 inline fn bitWidth(x: anytype) u16 {
14 const T = @TypeOf(x);
15 return switch (T) {
16 f80, f128, c_longdouble => @typeInfo(T).float.bits,
17 else => @compileError("Unsupported type: " ++ @typeName(T) ++ "\nPass a `f80`, `f128`, or `c_longdouble`."),
18 };
19 }
20
21 /// Returns the sign + exponent bits of a `long double`.
22 pub fn signExponent(x: anytype) u16 {
23 const bit_width = bitWidth(x);
24 switch (bit_width) {
25 80 => {
26 const bits: U80 = @bitCast(x);
27 return @intCast(bits >> 64);
28 },
29 128 => {
30 const bits: u128 = @bitCast(x);
31 return @intCast(bits >> 112);
32 },
33 // `c_longdouble` can have <80 bits on some targets, we want to error on that
34 else => @compileError(std.fmt.comptimePrint("`signExponent` supports floats of only `80` and `128` bit width, got bit width: {d}", .{bit_width})),
35 }
36 }
37
38 test "signExponent" {
39 try expectEqual(signExponent(@as(f80, -0.0)), 0x8000);
40 try expectEqual(signExponent(@as(f128, 0.0)), 0x0000);
41 try expectEqual(signExponent(@as(f128, 42.0)), 0x4004);
42 try expectEqual(signExponent(nan(c_longdouble)), 0x7FFF);
43 }
44
45 /// Takes the top 16 bits of a `long double`'s mantissa.
46 pub fn mantissaTop(x: anytype) u16 {
47 const bit_width = bitWidth(x);
48 switch (bit_width) {
49 80 => {
50 const bits: U80 = @bitCast(x);
51 return @intCast((bits >> 48) & 0xFFFF);
52 },
53 128 => {
54 const bits: u128 = @bitCast(x);
55 return @intCast((bits >> 96) & 0xFFFF);
56 },
57 // `c_longdouble` can have <80 bits on some targets, we want to error on that
58 else => @compileError(std.fmt.comptimePrint("`mantissaTop` supports floats of only `80` and `128` bit width, got bit width: {d}", .{bit_width})),
59 }
60 }
61
62 test "mantissaTop" {
63 try expectEqual(mantissaTop(@as(f80, -0.0)), 0x0000);
64 try expectEqual(mantissaTop(nan(f128)), 0x8000);
65 try expectEqual(mantissaTop(@as(f128, 42.0)), 0x5000);
66 }
67};
68
69pub fn FloatRepr(comptime Float: type) type {
70 const fractional_bits = floatFractionalBits(Float);
71 const exponent_bits = floatExponentBits(Float);
72 return packed struct {
73 const Repr = @This();
74
75 mantissa: StoredMantissa,
76 exponent: BiasedExponent,
77 sign: std.math.Sign,
78
79 pub const StoredMantissa = @Int(.unsigned, floatMantissaBits(Float));
80 pub const Mantissa = @Int(.unsigned, 1 + fractional_bits);
81 pub const Exponent = @Int(.signed, exponent_bits);
82 pub const BiasedExponent = enum(@Int(.unsigned, exponent_bits)) {
83 denormal = 0,
84 min_normal = 1,
85 zero = (1 << (exponent_bits - 1)) - 1,
86 max_normal = (1 << exponent_bits) - 2,
87 infinite = (1 << exponent_bits) - 1,
88 _,
89
90 pub const Int = @typeInfo(BiasedExponent).@"enum".tag_type;
91
92 pub fn unbias(biased: BiasedExponent) Exponent {
93 switch (biased) {
94 .denormal => unreachable,
95 else => return @bitCast(@backingInt(biased) -% @backingInt(BiasedExponent.zero)),
96 .infinite => unreachable,
97 }
98 }
99
100 pub fn bias(unbiased: Exponent) BiasedExponent {
101 return @fromBackingInt(@intCast(@backingInt(BiasedExponent.zero) +% @as(Int, @bitCast(unbiased))));
102 }
103 };
104
105 pub const Normalized = struct {
106 fraction: Fraction,
107 exponent: Normalized.Exponent,
108
109 pub const Fraction = @Int(.unsigned, fractional_bits);
110 pub const Exponent = @Int(.signed, 1 + exponent_bits);
111
112 /// This currently truncates denormal values, which needs to be fixed before this can be used to
113 /// produce a rounded value.
114 pub fn reconstruct(normalized: Normalized, sign: std.math.Sign) Float {
115 if (normalized.exponent > comptime BiasedExponent.max_normal.unbias()) return @bitCast(Repr{
116 .mantissa = 0,
117 .exponent = .infinite,
118 .sign = sign,
119 });
120 const mantissa = @as(Mantissa, 1 << fractional_bits) | normalized.fraction;
121 if (normalized.exponent < comptime BiasedExponent.min_normal.unbias()) return @bitCast(Repr{
122 .mantissa = @truncate(std.math.shr(
123 Mantissa,
124 mantissa,
125 (comptime BiasedExponent.min_normal.unbias()) - normalized.exponent,
126 )),
127 .exponent = .denormal,
128 .sign = sign,
129 });
130 return @bitCast(Repr{
131 .mantissa = @truncate(mantissa),
132 .exponent = .bias(@intCast(normalized.exponent)),
133 .sign = sign,
134 });
135 }
136 };
137
138 pub const Classified = union(enum) { normalized: Normalized, infinity, nan, invalid };
139 fn classify(repr: Repr) Classified {
140 return switch (repr.exponent) {
141 .denormal => {
142 const mantissa: Mantissa = repr.mantissa;
143 const shift = @clz(mantissa);
144 return .{ .normalized = .{
145 .fraction = @truncate(mantissa << shift),
146 .exponent = @as(Normalized.Exponent, comptime BiasedExponent.min_normal.unbias()) - shift,
147 } };
148 },
149 else => if (repr.mantissa <= std.math.maxInt(Normalized.Fraction)) .{ .normalized = .{
150 .fraction = @intCast(repr.mantissa),
151 .exponent = repr.exponent.unbias(),
152 } } else .invalid,
153 .infinite => switch (repr.mantissa) {
154 0 => .infinity,
155 else => .nan,
156 },
157 };
158 }
159 };
160}
161
162/// Creates a raw "1.0" mantissa for floating point type T. Used to dedupe f80 logic.
163inline fn mantissaOne(comptime T: type) comptime_int {
164 return if (@typeInfo(T).float.bits == 80) 1 << floatFractionalBits(T) else 0;
165}
166
167/// Creates floating point type T from an unbiased exponent and raw mantissa.
168inline fn reconstructFloat(comptime T: type, comptime exponent: comptime_int, comptime mantissa: comptime_int) T {
169 const TBits = @Int(.unsigned, @bitSizeOf(T));
170 const biased_exponent = @as(TBits, exponent + floatExponentMax(T));
171 return @as(T, @bitCast((biased_exponent << floatMantissaBits(T)) | @as(TBits, mantissa)));
172}
173
174/// Returns the number of bits in the exponent of floating point type T.
175pub inline fn floatExponentBits(comptime T: type) comptime_int {
176 comptime assert(@typeInfo(T) == .float);
177
178 return switch (@typeInfo(T).float.bits) {
179 16 => 5,
180 32 => 8,
181 64 => 11,
182 80 => 15,
183 128 => 15,
184 else => @compileError("unknown floating point type " ++ @typeName(T)),
185 };
186}
187
188/// Returns the number of bits in the mantissa of floating point type T.
189pub inline fn floatMantissaBits(comptime T: type) comptime_int {
190 comptime assert(@typeInfo(T) == .float);
191
192 return switch (@typeInfo(T).float.bits) {
193 16 => 10,
194 32 => 23,
195 64 => 52,
196 80 => 64,
197 128 => 112,
198 else => @compileError("unknown floating point type " ++ @typeName(T)),
199 };
200}
201
202/// Returns the number of fractional bits in the mantissa of floating point type T.
203pub inline fn floatFractionalBits(comptime T: type) comptime_int {
204 comptime assert(@typeInfo(T) == .float);
205
206 // standard IEEE floats have an implicit 0.m or 1.m integer part
207 // f80 is special and has an explicitly stored bit in the MSB
208 // this function corresponds to `MANT_DIG - 1' from C
209 return switch (@typeInfo(T).float.bits) {
210 16 => 10,
211 32 => 23,
212 64 => 52,
213 80 => 63,
214 128 => 112,
215 else => @compileError("unknown floating point type " ++ @typeName(T)),
216 };
217}
218
219/// Returns the minimum exponent that can represent
220/// a normalised value in floating point type T.
221pub inline fn floatExponentMin(comptime T: type) comptime_int {
222 return -floatExponentMax(T) + 1;
223}
224
225/// Returns the maximum exponent that can represent
226/// a normalised value in floating point type T.
227pub inline fn floatExponentMax(comptime T: type) comptime_int {
228 return (1 << (floatExponentBits(T) - 1)) - 1;
229}
230
231/// Returns the smallest subnormal number representable in floating point type T.
232pub inline fn floatTrueMin(comptime T: type) T {
233 return reconstructFloat(T, floatExponentMin(T) - 1, 1);
234}
235
236/// Returns the smallest normal number representable in floating point type T.
237pub inline fn floatMin(comptime T: type) T {
238 return reconstructFloat(T, floatExponentMin(T), mantissaOne(T));
239}
240
241/// Returns the largest normal number representable in floating point type T.
242pub inline fn floatMax(comptime T: type) T {
243 const all1s_mantissa = (1 << floatMantissaBits(T)) - 1;
244 return reconstructFloat(T, floatExponentMax(T), all1s_mantissa);
245}
246
247/// Returns the machine epsilon of floating point type T.
248pub inline fn floatEps(comptime T: type) T {
249 return reconstructFloat(T, -floatFractionalBits(T), mantissaOne(T));
250}
251
252/// Returns the local epsilon of floating point type T.
253pub inline fn floatEpsAt(comptime T: type, x: T) T {
254 switch (@typeInfo(T)) {
255 .float => |F| {
256 const U: type = @Int(.unsigned, F.bits);
257 const u: U = @bitCast(x);
258 const y: T = @bitCast(u ^ 1);
259 return @abs(x - y);
260 },
261 else => @compileError("floatEpsAt only supports floats"),
262 }
263}
264
265/// Returns the inf value for a floating point `Type`.
266pub inline fn inf(comptime Type: type) Type {
267 const RuntimeType = switch (Type) {
268 else => Type,
269 comptime_float => f128, // any float type will do
270 };
271 return reconstructFloat(RuntimeType, floatExponentMax(RuntimeType) + 1, mantissaOne(RuntimeType));
272}
273
274/// Returns the canonical quiet NaN representation for a floating point `Type`.
275pub inline fn nan(comptime Type: type) Type {
276 const RuntimeType = switch (Type) {
277 else => Type,
278 comptime_float => f128, // any float type will do
279 };
280 return reconstructFloat(
281 RuntimeType,
282 floatExponentMax(RuntimeType) + 1,
283 mantissaOne(RuntimeType) | 1 << (floatFractionalBits(RuntimeType) - 1),
284 );
285}
286
287/// Returns a signalling NaN representation for a floating point `Type`.
288///
289/// TODO: LLVM is known to miscompile on some architectures to quiet NaN -
290/// this is tracked by https://github.com/ziglang/zig/issues/14366
291pub inline fn snan(comptime Type: type) Type {
292 const RuntimeType = switch (Type) {
293 else => Type,
294 comptime_float => f128, // any float type will do
295 };
296 return reconstructFloat(
297 RuntimeType,
298 floatExponentMax(RuntimeType) + 1,
299 mantissaOne(RuntimeType) | 1 << (floatFractionalBits(RuntimeType) - 2),
300 );
301}
302
303fn floatBits(comptime Type: type) !void {
304 // (1 +) for the sign bit, since it is separate from the other bits
305 const size = 1 + floatExponentBits(Type) + floatMantissaBits(Type);
306 try expect(@bitSizeOf(Type) == size);
307 try expect(floatFractionalBits(Type) <= floatMantissaBits(Type));
308
309 // for machine epsilon, assert expmin <= -prec <= expmax
310 try expect(floatExponentMin(Type) <= -floatFractionalBits(Type));
311 try expect(-floatFractionalBits(Type) <= floatExponentMax(Type));
312}
313test floatBits {
314 try floatBits(f16);
315 try floatBits(f32);
316 try floatBits(f64);
317 try floatBits(f80);
318 try floatBits(f128);
319 try floatBits(c_longdouble);
320}
321
322test inf {
323 const inf_u16: u16 = 0x7C00;
324 const inf_u32: u32 = 0x7F800000;
325 const inf_u64: u64 = 0x7FF0000000000000;
326 const inf_u80: u80 = 0x7FFF8000000000000000;
327 const inf_u128: u128 = 0x7FFF0000000000000000000000000000;
328 try expectEqual(inf_u16, @as(u16, @bitCast(inf(f16))));
329 try expectEqual(inf_u32, @as(u32, @bitCast(inf(f32))));
330 try expectEqual(inf_u64, @as(u64, @bitCast(inf(f64))));
331 try expectEqual(inf_u80, @as(u80, @bitCast(inf(f80))));
332 try expectEqual(inf_u128, @as(u128, @bitCast(inf(f128))));
333}
334
335test nan {
336 const qnan_u16: u16 = 0x7E00;
337 const qnan_u32: u32 = 0x7FC00000;
338 const qnan_u64: u64 = 0x7FF8000000000000;
339 const qnan_u80: u80 = 0x7FFFC000000000000000;
340 const qnan_u128: u128 = 0x7FFF8000000000000000000000000000;
341 try expectEqual(qnan_u16, @as(u16, @bitCast(nan(f16))));
342 try expectEqual(qnan_u32, @as(u32, @bitCast(nan(f32))));
343 try expectEqual(qnan_u64, @as(u64, @bitCast(nan(f64))));
344 try expectEqual(qnan_u80, @as(u80, @bitCast(nan(f80))));
345 try expectEqual(qnan_u128, @as(u128, @bitCast(nan(f128))));
346}
347
348test snan {
349 const snan_u16: u16 = 0x7D00;
350 const snan_u32: u32 = 0x7FA00000;
351 const snan_u64: u64 = 0x7FF4000000000000;
352 const snan_u80: u80 = 0x7FFFA000000000000000;
353 const snan_u128: u128 = 0x7FFF4000000000000000000000000000;
354 try expectEqual(snan_u16, @as(u16, @bitCast(snan(f16))));
355 try expectEqual(snan_u32, @as(u32, @bitCast(snan(f32))));
356 try expectEqual(snan_u64, @as(u64, @bitCast(snan(f64))));
357 try expectEqual(snan_u80, @as(u80, @bitCast(snan(f80))));
358 try expectEqual(snan_u128, @as(u128, @bitCast(snan(f128))));
359}