authorgravatar for andrew@ziglang.orgAndrew Kelley <andrew@ziglang.org> 2022-04-19 13:04:31-04:00
committergravatar for noreply@github.comGitHub <noreply@github.com> 2022-04-19 13:04:31-04:00
logf3b3d7f20a51b85961e65257332842ad4030ece8
treee762437aeeb4bf5fd9195a25b84aa3553768a37b
parentedb4a07d4ddf01c7f4d589b0f427ba6b9e03134b
parentb2950866b18916fe5e6e458b62ec294244fa2c2f
signaturebadge-question-mark Signed by PGP key 4AEE18F83AFDEB23

Merge pull request #11468 from topolarity/f80-mul

compiler_rt: Implement softfloat multiply for `f80`

5 files changed, 264 insertions(+), 242 deletions(-)

lib/std/special/compiler_rt.zig+7-4
......@@ -226,23 +226,26 @@ comptime {
226226 @export(__addsf3, .{ .name = "__addsf3", .linkage = linkage });
227227 const __adddf3 = @import("compiler_rt/addXf3.zig").__adddf3;
228228 @export(__adddf3, .{ .name = "__adddf3", .linkage = linkage });
229 const __addtf3 = @import("compiler_rt/addXf3.zig").__addtf3;
230 @export(__addtf3, .{ .name = "__addtf3", .linkage = linkage });
231229 const __addxf3 = @import("compiler_rt/addXf3.zig").__addxf3;
232230 @export(__addxf3, .{ .name = "__addxf3", .linkage = linkage });
231 const __addtf3 = @import("compiler_rt/addXf3.zig").__addtf3;
232 @export(__addtf3, .{ .name = "__addtf3", .linkage = linkage });
233
233234 const __subsf3 = @import("compiler_rt/addXf3.zig").__subsf3;
234235 @export(__subsf3, .{ .name = "__subsf3", .linkage = linkage });
235236 const __subdf3 = @import("compiler_rt/addXf3.zig").__subdf3;
236237 @export(__subdf3, .{ .name = "__subdf3", .linkage = linkage });
237 const __subtf3 = @import("compiler_rt/addXf3.zig").__subtf3;
238 @export(__subtf3, .{ .name = "__subtf3", .linkage = linkage });
239238 const __subxf3 = @import("compiler_rt/addXf3.zig").__subxf3;
240239 @export(__subxf3, .{ .name = "__subxf3", .linkage = linkage });
240 const __subtf3 = @import("compiler_rt/addXf3.zig").__subtf3;
241 @export(__subtf3, .{ .name = "__subtf3", .linkage = linkage });
241242
242243 const __mulsf3 = @import("compiler_rt/mulXf3.zig").__mulsf3;
243244 @export(__mulsf3, .{ .name = "__mulsf3", .linkage = linkage });
244245 const __muldf3 = @import("compiler_rt/mulXf3.zig").__muldf3;
245246 @export(__muldf3, .{ .name = "__muldf3", .linkage = linkage });
247 const __mulxf3 = @import("compiler_rt/mulXf3.zig").__mulxf3;
248 @export(__mulxf3, .{ .name = "__mulxf3", .linkage = linkage });
246249 const __multf3 = @import("compiler_rt/mulXf3.zig").__multf3;
247250 @export(__multf3, .{ .name = "__multf3", .linkage = linkage });
248251
lib/std/special/compiler_rt/addXf3.zig+33-188
......@@ -3,6 +3,7 @@
33// https://github.com/llvm/llvm-project/blob/02d85149a05cb1f6dc49f0ba7a2ceca53718ae17/compiler-rt/lib/builtins/fp_add_impl.inc
44
55const std = @import("std");
6const math = std.math;
67const builtin = @import("builtin");
78const compiler_rt = @import("../compiler_rt.zig");
89
......@@ -14,6 +15,16 @@ pub fn __adddf3(a: f64, b: f64) callconv(.C) f64 {
1415 return addXf3(f64, a, b);
1516}
1617
18pub fn __addxf3(a: f80, b: f80) callconv(.C) f80 {
19 return addXf3(f80, a, b);
20}
21
22pub fn __subxf3(a: f80, b: f80) callconv(.C) f80 {
23 var b_rep = std.math.break_f80(b);
24 b_rep.exp ^= 0x8000;
25 return __addxf3(a, std.math.make_f80(b_rep));
26}
27
1728pub fn __addtf3(a: f128, b: f128) callconv(.C) f128 {
1829 return addXf3(f128, a, b);
1930}
......@@ -58,10 +69,10 @@ fn normalize(comptime T: type, significand: *std.meta.Int(.unsigned, @typeInfo(T
5869 const bits = @typeInfo(T).Float.bits;
5970 const Z = std.meta.Int(.unsigned, bits);
6071 const S = std.meta.Int(.unsigned, bits - @clz(Z, @as(Z, bits) - 1));
61 const significandBits = std.math.floatMantissaBits(T);
62 const implicitBit = @as(Z, 1) << significandBits;
72 const fractionalBits = math.floatFractionalBits(T);
73 const integerBit = @as(Z, 1) << fractionalBits;
6374
64 const shift = @clz(std.meta.Int(.unsigned, bits), significand.*) - @clz(Z, implicitBit);
75 const shift = @clz(std.meta.Int(.unsigned, bits), significand.*) - @clz(Z, integerBit);
6576 significand.* <<= @intCast(S, shift);
6677 return 1 - shift;
6778}
......@@ -73,26 +84,26 @@ fn addXf3(comptime T: type, a: T, b: T) T {
7384 const S = std.meta.Int(.unsigned, bits - @clz(Z, @as(Z, bits) - 1));
7485
7586 const typeWidth = bits;
76 const significandBits = std.math.floatMantissaBits(T);
77 const exponentBits = std.math.floatExponentBits(T);
87 const significandBits = math.floatMantissaBits(T);
88 const fractionalBits = math.floatFractionalBits(T);
89 const exponentBits = math.floatExponentBits(T);
7890
7991 const signBit = (@as(Z, 1) << (significandBits + exponentBits));
8092 const maxExponent = ((1 << exponentBits) - 1);
8193
82 const implicitBit = (@as(Z, 1) << significandBits);
83 const quietBit = implicitBit >> 1;
84 const significandMask = implicitBit - 1;
94 const integerBit = (@as(Z, 1) << fractionalBits);
95 const quietBit = integerBit >> 1;
96 const significandMask = (@as(Z, 1) << significandBits) - 1;
8597
8698 const absMask = signBit - 1;
87 const exponentMask = absMask ^ significandMask;
88 const qnanRep = exponentMask | quietBit;
99 const qnanRep = @bitCast(Z, math.nan(T)) | quietBit;
89100
90101 var aRep = @bitCast(Z, a);
91102 var bRep = @bitCast(Z, b);
92103 const aAbs = aRep & absMask;
93104 const bAbs = bRep & absMask;
94105
95 const infRep = @bitCast(Z, std.math.inf(T));
106 const infRep = @bitCast(Z, math.inf(T));
96107
97108 // Detect if a or b is zero, infinity, or NaN.
98109 if (aAbs -% @as(Z, 1) >= infRep - @as(Z, 1) or
......@@ -157,12 +168,12 @@ fn addXf3(comptime T: type, a: T, b: T) T {
157168 // implicit significand bit. (If we fell through from the denormal path it
158169 // was already set by normalize( ), but setting it twice won't hurt
159170 // anything.)
160 aSignificand = (aSignificand | implicitBit) << 3;
161 bSignificand = (bSignificand | implicitBit) << 3;
171 aSignificand = (aSignificand | integerBit) << 3;
172 bSignificand = (bSignificand | integerBit) << 3;
162173
163174 // Shift the significand of b by the difference in exponents, with a sticky
164175 // bottom bit to get rounding correct.
165 const @"align" = @intCast(Z, aExponent - bExponent);
176 const @"align" = @intCast(u32, aExponent - bExponent);
166177 if (@"align" != 0) {
167178 if (@"align" < typeWidth) {
168179 const sticky = if (bSignificand << @intCast(S, typeWidth - @"align") != 0) @as(Z, 1) else 0;
......@@ -178,8 +189,8 @@ fn addXf3(comptime T: type, a: T, b: T) T {
178189
179190 // If partial cancellation occured, we need to left-shift the result
180191 // and adjust the exponent:
181 if (aSignificand < implicitBit << 3) {
182 const shift = @intCast(i32, @clz(Z, aSignificand)) - @intCast(i32, @clz(std.meta.Int(.unsigned, bits), implicitBit << 3));
192 if (aSignificand < integerBit << 3) {
193 const shift = @intCast(i32, @clz(Z, aSignificand)) - @intCast(i32, @clz(std.meta.Int(.unsigned, bits), integerBit << 3));
183194 aSignificand <<= @intCast(S, shift);
184195 aExponent -= shift;
185196 }
......@@ -188,7 +199,7 @@ fn addXf3(comptime T: type, a: T, b: T) T {
188199
189200 // If the addition carried up, we need to right-shift the result and
190201 // adjust the exponent:
191 if (aSignificand & (implicitBit << 4) != 0) {
202 if (aSignificand & (integerBit << 4) != 0) {
192203 const sticky = aSignificand & 1;
193204 aSignificand = aSignificand >> 1 | sticky;
194205 aExponent += 1;
......@@ -210,7 +221,7 @@ fn addXf3(comptime T: type, a: T, b: T) T {
210221 // Low three bits are round, guard, and sticky.
211222 const roundGuardSticky = aSignificand & 0x7;
212223
213 // Shift the significand into place, and mask off the implicit bit.
224 // Shift the significand into place, and mask off the integer bit, if it's implicit.
214225 var result = (aSignificand >> 3) & significandMask;
215226
216227 // Insert the exponent and sign.
......@@ -222,178 +233,12 @@ fn addXf3(comptime T: type, a: T, b: T) T {
222233 if (roundGuardSticky > 0x4) result += 1;
223234 if (roundGuardSticky == 0x4) result += result & 1;
224235
225 return @bitCast(T, result);
226}
227
228fn normalize_f80(exp: *i32, significand: *u80) void {
229 const shift = @clz(u64, @truncate(u64, significand.*));
230 significand.* = (significand.* << shift);
231 exp.* += -@as(i8, shift);
232}
233
234pub fn __addxf3(a: f80, b: f80) callconv(.C) f80 {
235 var a_rep = std.math.break_f80(a);
236 var b_rep = std.math.break_f80(b);
237 var a_exp: i32 = a_rep.exp & 0x7FFF;
238 var b_exp: i32 = b_rep.exp & 0x7FFF;
239
240 const significand_bits = std.math.floatMantissaBits(f80);
241 const int_bit = 0x8000000000000000;
242 const significand_mask = 0x7FFFFFFFFFFFFFFF;
243 const qnan_bit = 0xC000000000000000;
244 const max_exp = 0x7FFF;
245 const sign_bit = 0x8000;
246
247 // Detect if a or b is infinity, or NaN.
248 if (a_exp == max_exp) {
249 if (a_rep.fraction ^ int_bit == 0) {
250 if (b_exp == max_exp and (b_rep.fraction ^ int_bit == 0)) {
251 // +/-infinity + -/+infinity = qNaN
252 return std.math.qnan_f80;
253 }
254 // +/-infinity + anything = +/- infinity
255 return a;
256 } else {
257 std.debug.assert(a_rep.fraction & significand_mask != 0);
258 // NaN + anything = qNaN
259 a_rep.fraction |= qnan_bit;
260 return std.math.make_f80(a_rep);
261 }
262 }
263 if (b_exp == max_exp) {
264 if (b_rep.fraction ^ int_bit == 0) {
265 // anything + +/-infinity = +/-infinity
266 return b;
267 } else {
268 std.debug.assert(b_rep.fraction & significand_mask != 0);
269 // anything + NaN = qNaN
270 b_rep.fraction |= qnan_bit;
271 return std.math.make_f80(b_rep);
272 }
273 }
274
275 const a_zero = (a_rep.fraction | @bitCast(u32, a_exp)) == 0;
276 const b_zero = (b_rep.fraction | @bitCast(u32, b_exp)) == 0;
277 if (a_zero) {
278 // zero + anything = anything
279 if (b_zero) {
280 // but we need to get the sign right for zero + zero
281 a_rep.exp &= b_rep.exp;
282 return std.math.make_f80(a_rep);
283 } else {
284 return b;
285 }
286 } else if (b_zero) {
287 // anything + zero = anything
288 return a;
289 }
290
291 var a_int: u80 = a_rep.fraction | (@as(u80, a_rep.exp & max_exp) << significand_bits);
292 var b_int: u80 = b_rep.fraction | (@as(u80, b_rep.exp & max_exp) << significand_bits);
293
294 // Swap a and b if necessary so that a has the larger absolute value.
295 if (b_int > a_int) {
296 const temp = a_rep;
297 a_rep = b_rep;
298 b_rep = temp;
299 }
300
301 // Extract the exponent and significand from the (possibly swapped) a and b.
302 a_exp = a_rep.exp & max_exp;
303 b_exp = b_rep.exp & max_exp;
304 a_int = a_rep.fraction;
305 b_int = b_rep.fraction;
306
307 // Normalize any denormals, and adjust the exponent accordingly.
308 normalize_f80(&a_exp, &a_int);
309 normalize_f80(&b_exp, &b_int);
310
311 // The sign of the result is the sign of the larger operand, a. If they
312 // have opposite signs, we are performing a subtraction; otherwise addition.
313 const result_sign = a_rep.exp & sign_bit;
314 const subtraction = (a_rep.exp ^ b_rep.exp) & sign_bit != 0;
315
316 // Shift the significands to give us round, guard and sticky, and or in the
317 // implicit significand bit. (If we fell through from the denormal path it
318 // was already set by normalize( ), but setting it twice won't hurt
319 // anything.)
320 a_int = a_int << 3;
321 b_int = b_int << 3;
322
323 // Shift the significand of b by the difference in exponents, with a sticky
324 // bottom bit to get rounding correct.
325 const @"align" = @intCast(u80, a_exp - b_exp);
326 if (@"align" != 0) {
327 if (@"align" < 80) {
328 const sticky = if (b_int << @intCast(u7, 80 - @"align") != 0) @as(u80, 1) else 0;
329 b_int = (b_int >> @truncate(u7, @"align")) | sticky;
330 } else {
331 b_int = 1; // sticky; b is known to be non-zero.
332 }
333 }
334 if (subtraction) {
335 a_int -= b_int;
336 // If a == -b, return +zero.
337 if (a_int == 0) return 0.0;
338
339 // If partial cancellation occurred, we need to left-shift the result
340 // and adjust the exponent:
341 if (a_int < int_bit << 3) {
342 const shift = @intCast(i32, @clz(u80, a_int)) - @intCast(i32, @clz(u80, @as(u80, int_bit) << 3));
343 a_int <<= @intCast(u7, shift);
344 a_exp -= shift;
345 }
346 } else { // addition
347 a_int += b_int;
348
349 // If the addition carried up, we need to right-shift the result and
350 // adjust the exponent:
351 if (a_int & (int_bit << 4) != 0) {
352 const sticky = a_int & 1;
353 a_int = a_int >> 1 | sticky;
354 a_exp += 1;
355 }
236 // Restore any explicit integer bit, if it was rounded off
237 if (significandBits != fractionalBits) {
238 if ((result >> significandBits) != 0) result |= integerBit;
356239 }
357240
358 // If we have overflowed the type, return +/- infinity:
359 if (a_exp >= max_exp) {
360 a_rep.exp = max_exp | result_sign;
361 a_rep.fraction = int_bit; // integer bit is set for +/-inf
362 return std.math.make_f80(a_rep);
363 }
364
365 if (a_exp <= 0) {
366 // Result is denormal before rounding; the exponent is zero and we
367 // need to shift the significand.
368 const shift = @intCast(u80, 1 - a_exp);
369 const sticky = if (a_int << @intCast(u7, 80 - shift) != 0) @as(u1, 1) else 0;
370 a_int = a_int >> @intCast(u7, shift | sticky);
371 a_exp = 0;
372 }
373
374 // Low three bits are round, guard, and sticky.
375 const round_guard_sticky = @truncate(u3, a_int);
376
377 // Shift the significand into place.
378 a_int = @truncate(u64, a_int >> 3);
379
380 // // Insert the exponent and sign.
381 a_int |= (@intCast(u80, a_exp) | result_sign) << significand_bits;
382
383 // Final rounding. The result may overflow to infinity, but that is the
384 // correct result in that case.
385 if (round_guard_sticky > 0x4) a_int += 1;
386 if (round_guard_sticky == 0x4) a_int += a_int & 1;
387
388 a_rep.fraction = @truncate(u64, a_int);
389 a_rep.exp = @truncate(u16, a_int >> significand_bits);
390 return std.math.make_f80(a_rep);
391}
392
393pub fn __subxf3(a: f80, b: f80) callconv(.C) f80 {
394 var b_rep = std.math.break_f80(b);
395 b_rep.exp ^= 0x8000;
396 return __addxf3(a, std.math.make_f80(b_rep));
241 return @bitCast(T, result);
397242}
398243
399244test {
lib/std/special/compiler_rt/addXf3_test.zig+74-4
......@@ -3,8 +3,9 @@
33// https://github.com/llvm/llvm-project/blob/02d85149a05cb1f6dc49f0ba7a2ceca53718ae17/compiler-rt/test/builtins/Unit/addtf3_test.c
44// https://github.com/llvm/llvm-project/blob/02d85149a05cb1f6dc49f0ba7a2ceca53718ae17/compiler-rt/test/builtins/Unit/subtf3_test.c
55
6const std = @import("std");
7const math = std.math;
68const qnan128 = @bitCast(f128, @as(u128, 0x7fff800000000000) << 64);
7const inf128 = @bitCast(f128, @as(u128, 0x7fff000000000000) << 64);
89
910const __addtf3 = @import("addXf3.zig").__addtf3;
1011
......@@ -37,13 +38,14 @@ test "addtf3" {
3738 try test__addtf3(@bitCast(f128, (@as(u128, 0x7fff000000000000) << 64) | @as(u128, 0x800030000000)), 0x1.23456789abcdefp+5, 0x7fff800000000000, 0x0);
3839
3940 // inf + inf = inf
40 try test__addtf3(inf128, inf128, 0x7fff000000000000, 0x0);
41 try test__addtf3(math.inf(f128), math.inf(f128), 0x7fff000000000000, 0x0);
4142
4243 // inf + any = inf
43 try test__addtf3(inf128, 0x1.2335653452436234723489432abcdefp+5, 0x7fff000000000000, 0x0);
44 try test__addtf3(math.inf(f128), 0x1.2335653452436234723489432abcdefp+5, 0x7fff000000000000, 0x0);
4445
4546 // any + any
4647 try test__addtf3(0x1.23456734245345543849abcdefp+5, 0x1.edcba52449872455634654321fp-1, 0x40042afc95c8b579, 0x61e58dd6c51eb77c);
48 try test__addtf3(0x1.edcba52449872455634654321fp-1, 0x1.23456734245345543849abcdefp+5, 0x40042afc95c8b579, 0x61e58dd6c51eb77c);
4749}
4850
4951const __subtf3 = @import("addXf3.zig").__subtf3;
......@@ -78,8 +80,76 @@ test "subtf3" {
7880 try test__subtf3(@bitCast(f128, (@as(u128, 0x7fff000000000000) << 64) | @as(u128, 0x800030000000)), 0x1.23456789abcdefp+5, 0x7fff800000000000, 0x0);
7981
8082 // inf - any = inf
81 try test__subtf3(inf128, 0x1.23456789abcdefp+5, 0x7fff000000000000, 0x0);
83 try test__subtf3(math.inf(f128), 0x1.23456789abcdefp+5, 0x7fff000000000000, 0x0);
8284
8385 // any + any
8486 try test__subtf3(0x1.234567829a3bcdef5678ade36734p+5, 0x1.ee9d7c52354a6936ab8d7654321fp-1, 0x40041b8af1915166, 0xa44a7bca780a166c);
87 try test__subtf3(0x1.ee9d7c52354a6936ab8d7654321fp-1, 0x1.234567829a3bcdef5678ade36734p+5, 0xc0041b8af1915166, 0xa44a7bca780a166c);
88}
89
90const __addxf3 = @import("addXf3.zig").__addxf3;
91const qnan80 = @bitCast(f80, @bitCast(u80, math.nan(f80)) | (1 << (math.floatFractionalBits(f80) - 1)));
92
93fn test__addxf3(a: f80, b: f80, expected: u80) !void {
94 const x = __addxf3(a, b);
95 const rep = @bitCast(u80, x);
96
97 if (rep == expected)
98 return;
99
100 if (math.isNan(@bitCast(f80, expected)) and math.isNan(x))
101 return; // We don't currently test NaN payload propagation
102
103 return error.TestFailed;
104}
105
106test "addxf3" {
107 // NaN + any = NaN
108 try test__addxf3(qnan80, 0x1.23456789abcdefp+5, @bitCast(u80, qnan80));
109 try test__addxf3(@bitCast(f80, @as(u80, 0x7fff_8000_8000_3000_0000)), 0x1.23456789abcdefp+5, @bitCast(u80, qnan80));
110
111 // any + NaN = NaN
112 try test__addxf3(0x1.23456789abcdefp+5, qnan80, @bitCast(u80, qnan80));
113 try test__addxf3(0x1.23456789abcdefp+5, @bitCast(f80, @as(u80, 0x7fff_8000_8000_3000_0000)), @bitCast(u80, qnan80));
114
115 // NaN + inf = NaN
116 try test__addxf3(qnan80, math.inf(f80), @bitCast(u80, qnan80));
117
118 // inf + NaN = NaN
119 try test__addxf3(math.inf(f80), qnan80, @bitCast(u80, qnan80));
120
121 // inf + inf = inf
122 try test__addxf3(math.inf(f80), math.inf(f80), @bitCast(u80, math.inf(f80)));
123
124 // inf + -inf = NaN
125 try test__addxf3(math.inf(f80), -math.inf(f80), @bitCast(u80, qnan80));
126
127 // -inf + inf = NaN
128 try test__addxf3(-math.inf(f80), math.inf(f80), @bitCast(u80, qnan80));
129
130 // inf + any = inf
131 try test__addxf3(math.inf(f80), 0x1.2335653452436234723489432abcdefp+5, @bitCast(u80, math.inf(f80)));
132
133 // any + inf = inf
134 try test__addxf3(0x1.2335653452436234723489432abcdefp+5, math.inf(f80), @bitCast(u80, math.inf(f80)));
135
136 // any + any
137 try test__addxf3(0x1.23456789abcdp+5, 0x1.dcba987654321p+5, 0x4005_BFFFFFFFFFFFC400);
138 try test__addxf3(0x1.23456734245345543849abcdefp+5, 0x1.edcba52449872455634654321fp-1, 0x4004_957E_4AE4_5ABC_B0F3);
139 try test__addxf3(0x1.ffff_ffff_ffff_fffcp+0, 0x1.0p-63, 0x3FFF_FFFFFFFFFFFFFFFF); // exact
140 try test__addxf3(0x1.ffff_ffff_ffff_fffep+0, 0x0.0p0, 0x3FFF_FFFFFFFFFFFFFFFF); // exact
141 try test__addxf3(0x1.ffff_ffff_ffff_fffcp+0, 0x1.4p-63, 0x3FFF_FFFFFFFFFFFFFFFF); // round down
142 try test__addxf3(0x1.ffff_ffff_ffff_fffcp+0, 0x1.8p-63, 0x4000_8000000000000000); // round up to even
143 try test__addxf3(0x1.ffff_ffff_ffff_fffcp+0, 0x1.cp-63, 0x4000_8000000000000000); // round up
144 try test__addxf3(0x1.ffff_ffff_ffff_fffcp+0, 0x2.0p-63, 0x4000_8000000000000000); // exact
145 try test__addxf3(0x1.ffff_ffff_ffff_fffcp+0, 0x2.1p-63, 0x4000_8000000000000000); // round down
146 try test__addxf3(0x1.ffff_ffff_ffff_fffcp+0, 0x3.0p-63, 0x4000_8000000000000000); // round down to even
147 try test__addxf3(0x1.ffff_ffff_ffff_fffcp+0, 0x3.1p-63, 0x4000_8000000000000001); // round up
148 try test__addxf3(0x1.ffff_ffff_ffff_fffcp+0, 0x4.0p-63, 0x4000_8000000000000001); // exact
149
150 try test__addxf3(0x1.0fff_ffff_ffff_fffep+0, 0x1.0p-63, 0x3FFF_8800000000000000); // exact
151 try test__addxf3(0x1.0fff_ffff_ffff_fffep+0, 0x1.7p-63, 0x3FFF_8800000000000000); // round down
152 try test__addxf3(0x1.0fff_ffff_ffff_fffep+0, 0x1.8p-63, 0x3FFF_8800000000000000); // round down to even
153 try test__addxf3(0x1.0fff_ffff_ffff_fffep+0, 0x1.9p-63, 0x3FFF_8800000000000001); // round up
154 try test__addxf3(0x1.0fff_ffff_ffff_fffep+0, 0x2.0p-63, 0x3FFF_8800000000000001); // exact
85155}
lib/std/special/compiler_rt/mulXf3.zig+83-46
......@@ -3,12 +3,16 @@
33// https://github.com/llvm/llvm-project/blob/2ffb1b0413efa9a24eb3c49e710e36f92e2cb50b/compiler-rt/lib/builtins/fp_mul_impl.inc
44
55const std = @import("std");
6const math = std.math;
67const builtin = @import("builtin");
78const compiler_rt = @import("../compiler_rt.zig");
89
910pub fn __multf3(a: f128, b: f128) callconv(.C) f128 {
1011 return mulXf3(f128, a, b);
1112}
13pub fn __mulxf3(a: f80, b: f80) callconv(.C) f80 {
14 return mulXf3(f80, a, b);
15}
1216pub fn __muldf3(a: f64, b: f64) callconv(.C) f64 {
1317 return mulXf3(f64, a, b);
1418}
......@@ -29,30 +33,36 @@ pub fn __aeabi_dmul(a: f64, b: f64) callconv(.C) f64 {
2933fn mulXf3(comptime T: type, a: T, b: T) T {
3034 @setRuntimeSafety(builtin.is_test);
3135 const typeWidth = @typeInfo(T).Float.bits;
36 const significandBits = math.floatMantissaBits(T);
37 const fractionalBits = math.floatFractionalBits(T);
38 const exponentBits = math.floatExponentBits(T);
39
3240 const Z = std.meta.Int(.unsigned, typeWidth);
3341
34 const significandBits = std.math.floatMantissaBits(T);
35 const exponentBits = std.math.floatExponentBits(T);
42 // ZSignificand is large enough to contain the significand, including an explicit integer bit
43 const ZSignificand = PowerOfTwoSignificandZ(T);
44 const ZSignificandBits = @typeInfo(ZSignificand).Int.bits;
3645
46 const roundBit = (1 << (ZSignificandBits - 1));
3747 const signBit = (@as(Z, 1) << (significandBits + exponentBits));
3848 const maxExponent = ((1 << exponentBits) - 1);
3949 const exponentBias = (maxExponent >> 1);
4050
41 const implicitBit = (@as(Z, 1) << significandBits);
42 const quietBit = implicitBit >> 1;
43 const significandMask = implicitBit - 1;
51 const integerBit = (@as(ZSignificand, 1) << fractionalBits);
52 const quietBit = integerBit >> 1;
53 const significandMask = (@as(Z, 1) << significandBits) - 1;
4454
4555 const absMask = signBit - 1;
46 const exponentMask = absMask ^ significandMask;
47 const qnanRep = exponentMask | quietBit;
48 const infRep = @bitCast(Z, std.math.inf(T));
56 const qnanRep = @bitCast(Z, math.nan(T)) | quietBit;
57 const infRep = @bitCast(Z, math.inf(T));
58 const minNormalRep = @bitCast(Z, math.floatMin(T));
4959
5060 const aExponent = @truncate(u32, (@bitCast(Z, a) >> significandBits) & maxExponent);
5161 const bExponent = @truncate(u32, (@bitCast(Z, b) >> significandBits) & maxExponent);
5262 const productSign: Z = (@bitCast(Z, a) ^ @bitCast(Z, b)) & signBit;
5363
54 var aSignificand: Z = @bitCast(Z, a) & significandMask;
55 var bSignificand: Z = @bitCast(Z, b) & significandMask;
64 var aSignificand: ZSignificand = @intCast(ZSignificand, @bitCast(Z, a) & significandMask);
65 var bSignificand: ZSignificand = @intCast(ZSignificand, @bitCast(Z, b) & significandMask);
5666 var scale: i32 = 0;
5767
5868 // Detect if a or b is zero, denormal, infinity, or NaN.
......@@ -93,38 +103,40 @@ fn mulXf3(comptime T: type, a: T, b: T) T {
93103 // one or both of a or b is denormal, the other (if applicable) is a
94104 // normal number. Renormalize one or both of a and b, and set scale to
95105 // include the necessary exponent adjustment.
96 if (aAbs < implicitBit) scale += normalize(T, &aSignificand);
97 if (bAbs < implicitBit) scale += normalize(T, &bSignificand);
106 if (aAbs < minNormalRep) scale += normalize(T, &aSignificand);
107 if (bAbs < minNormalRep) scale += normalize(T, &bSignificand);
98108 }
99109
100110 // Or in the implicit significand bit. (If we fell through from the
101111 // denormal path it was already set by normalize( ), but setting it twice
102112 // won't hurt anything.)
103 aSignificand |= implicitBit;
104 bSignificand |= implicitBit;
113 aSignificand |= integerBit;
114 bSignificand |= integerBit;
105115
106116 // Get the significand of a*b. Before multiplying the significands, shift
107117 // one of them left to left-align it in the field. Thus, the product will
108118 // have (exponentBits + 2) integral digits, all but two of which must be
109119 // zero. Normalizing this result is just a conditional left-shift by one
110120 // and bumping the exponent accordingly.
111 var productHi: Z = undefined;
112 var productLo: Z = undefined;
113 wideMultiply(Z, aSignificand, bSignificand << exponentBits, &productHi, &productLo);
121 var productHi: ZSignificand = undefined;
122 var productLo: ZSignificand = undefined;
123 const left_align_shift = ZSignificandBits - fractionalBits - 1;
124 wideMultiply(ZSignificand, aSignificand, bSignificand << left_align_shift, &productHi, &productLo);
114125
115 var productExponent: i32 = @bitCast(i32, aExponent +% bExponent) -% exponentBias +% scale;
126 var productExponent: i32 = @intCast(i32, aExponent + bExponent) - exponentBias + scale;
116127
117128 // Normalize the significand, adjust exponent if needed.
118 if ((productHi & implicitBit) != 0) {
129 if ((productHi & integerBit) != 0) {
119130 productExponent +%= 1;
120131 } else {
121 productHi = (productHi << 1) | (productLo >> (typeWidth - 1));
132 productHi = (productHi << 1) | (productLo >> (ZSignificandBits - 1));
122133 productLo = productLo << 1;
123134 }
124135
125136 // If we have overflowed the type, return +/- infinity.
126137 if (productExponent >= maxExponent) return @bitCast(T, infRep | productSign);
127138
139 var result: Z = undefined;
128140 if (productExponent <= 0) {
129141 // Result is denormal before rounding
130142 //
......@@ -133,35 +145,49 @@ fn mulXf3(comptime T: type, a: T, b: T) T {
133145 // handle this case separately, but we make it a special case to
134146 // simplify the shift logic.
135147 const shift: u32 = @truncate(u32, @as(Z, 1) -% @bitCast(u32, productExponent));
136 if (shift >= typeWidth) return @bitCast(T, productSign);
148 if (shift >= ZSignificandBits) return @bitCast(T, productSign);
137149
138150 // Otherwise, shift the significand of the result so that the round
139151 // bit is the high bit of productLo.
140 wideRightShiftWithSticky(Z, &productHi, &productLo, shift);
152 const sticky = wideShrWithTruncation(ZSignificand, &productHi, &productLo, shift);
153 productLo |= @boolToInt(sticky);
154 result = productHi;
141155 } else {
142156 // Result is normal before rounding; insert the exponent.
143 productHi &= significandMask;
144 productHi |= @as(Z, @bitCast(u32, productExponent)) << significandBits;
157 result = productHi & significandMask;
158 result |= @intCast(Z, productExponent) << significandBits;
145159 }
146160
147 // Insert the sign of the result:
148 productHi |= productSign;
149
150161 // Final rounding. The final result may overflow to infinity, or underflow
151162 // to zero, but those are the correct results in those cases. We use the
152163 // default IEEE-754 round-to-nearest, ties-to-even rounding mode.
153 if (productLo > signBit) productHi +%= 1;
154 if (productLo == signBit) productHi +%= productHi & 1;
155 return @bitCast(T, productHi);
164 if (productLo > roundBit) result +%= 1;
165 if (productLo == roundBit) result +%= result & 1;
166
167 // Restore any explicit integer bit, if it was rounded off
168 if (significandBits != fractionalBits) {
169 if ((result >> significandBits) != 0) result |= integerBit;
170 }
171
172 // Insert the sign of the result:
173 result |= productSign;
174
175 return @bitCast(T, result);
156176}
157177
158178fn wideMultiply(comptime Z: type, a: Z, b: Z, hi: *Z, lo: *Z) void {
159179 @setRuntimeSafety(builtin.is_test);
160180 switch (Z) {
181 u16 => {
182 // 16x16 --> 32 bit multiply
183 const product = @as(u32, a) * @as(u32, b);
184 hi.* = @intCast(u16, product >> 16);
185 lo.* = @truncate(u16, product);
186 },
161187 u32 => {
162188 // 32x32 --> 64 bit multiply
163189 const product = @as(u64, a) * @as(u64, b);
164 hi.* = @truncate(u32, product >> 32);
190 hi.* = @intCast(u32, product >> 32);
165191 lo.* = @truncate(u32, product);
166192 },
167193 u64 => {
......@@ -170,7 +196,7 @@ fn wideMultiply(comptime Z: type, a: Z, b: Z, hi: *Z, lo: *Z) void {
170196 return @truncate(u32, x);
171197 }
172198 fn hiWord(x: u64) u64 {
173 return @truncate(u32, x >> 32);
199 return @intCast(u32, x >> 32);
174200 }
175201 };
176202 // 64x64 -> 128 wide multiply for platforms that don't have such an operation;
......@@ -264,34 +290,45 @@ fn wideMultiply(comptime Z: type, a: Z, b: Z, hi: *Z, lo: *Z) void {
264290 }
265291}
266292
267fn normalize(comptime T: type, significand: *std.meta.Int(.unsigned, @typeInfo(T).Float.bits)) i32 {
293/// Returns a power-of-two integer type that is large enough to contain
294/// the significand of T, including an explicit integer bit
295fn PowerOfTwoSignificandZ(comptime T: type) type {
296 const bits = math.ceilPowerOfTwoAssert(u16, math.floatFractionalBits(T) + 1);
297 return std.meta.Int(.unsigned, bits);
298}
299
300fn normalize(comptime T: type, significand: *PowerOfTwoSignificandZ(T)) i32 {
268301 @setRuntimeSafety(builtin.is_test);
269 const Z = std.meta.Int(.unsigned, @typeInfo(T).Float.bits);
270 const significandBits = std.math.floatMantissaBits(T);
271 const implicitBit = @as(Z, 1) << significandBits;
302 const Z = PowerOfTwoSignificandZ(T);
303 const integerBit = @as(Z, 1) << math.floatFractionalBits(T);
272304
273 const shift = @clz(Z, significand.*) - @clz(Z, implicitBit);
274 significand.* <<= @intCast(std.math.Log2Int(Z), shift);
305 const shift = @clz(Z, significand.*) - @clz(Z, integerBit);
306 significand.* <<= @intCast(math.Log2Int(Z), shift);
275307 return @as(i32, 1) - shift;
276308}
277309
278fn wideRightShiftWithSticky(comptime Z: type, hi: *Z, lo: *Z, count: u32) void {
310// Returns `true` if the right shift is inexact (i.e. any bit shifted out is non-zero)
311//
312// This is analogous to an shr version of `@shlWithOverflow`
313fn wideShrWithTruncation(comptime Z: type, hi: *Z, lo: *Z, count: u32) bool {
279314 @setRuntimeSafety(builtin.is_test);
280315 const typeWidth = @typeInfo(Z).Int.bits;
281 const S = std.math.Log2Int(Z);
316 const S = math.Log2Int(Z);
317 var inexact = false;
282318 if (count < typeWidth) {
283 const sticky = @boolToInt((lo.* << @intCast(S, typeWidth -% count)) != 0);
284 lo.* = (hi.* << @intCast(S, typeWidth -% count)) | (lo.* >> @intCast(S, count)) | sticky;
319 inexact = (lo.* << @intCast(S, typeWidth -% count)) != 0;
320 lo.* = (hi.* << @intCast(S, typeWidth -% count)) | (lo.* >> @intCast(S, count));
285321 hi.* = hi.* >> @intCast(S, count);
286322 } else if (count < 2 * typeWidth) {
287 const sticky = @boolToInt((hi.* << @intCast(S, 2 * typeWidth -% count) | lo.*) != 0);
288 lo.* = hi.* >> @intCast(S, count -% typeWidth) | sticky;
323 inexact = (hi.* << @intCast(S, 2 * typeWidth -% count) | lo.*) != 0;
324 lo.* = hi.* >> @intCast(S, count -% typeWidth);
289325 hi.* = 0;
290326 } else {
291 const sticky = @boolToInt((hi.* | lo.*) != 0);
292 lo.* = sticky;
327 inexact = (hi.* | lo.*) != 0;
328 lo.* = 0;
293329 hi.* = 0;
294330 }
331 return inexact;
295332}
296333
297334test {
lib/std/special/compiler_rt/mulXf3_test.zig+67
......@@ -2,10 +2,15 @@
22//
33// https://github.com/llvm/llvm-project/blob/2ffb1b0413efa9a24eb3c49e710e36f92e2cb50b/compiler-rt/test/builtins/Unit/multf3_test.c
44
5const std = @import("std");
6const math = std.math;
57const qnan128 = @bitCast(f128, @as(u128, 0x7fff800000000000) << 64);
68const inf128 = @bitCast(f128, @as(u128, 0x7fff000000000000) << 64);
79
810const __multf3 = @import("mulXf3.zig").__multf3;
11const __mulxf3 = @import("mulXf3.zig").__mulxf3;
12const __muldf3 = @import("mulXf3.zig").__muldf3;
13const __mulsf3 = @import("mulXf3.zig").__mulsf3;
914
1015// return true if equal
1116// use two 64-bit integers intead of one 128-bit integer
......@@ -97,4 +102,66 @@ test "multf3" {
97102 0x3f90000000000000,
98103 0x0,
99104 );
105
106 try test__multf3(0x1.0000_0000_0000_0000_0000_0000_0001p+0, 0x1.8p+5, 0x4004_8000_0000_0000, 0x0000_0000_0000_0002);
107 try test__multf3(0x1.0000_0000_0000_0000_0000_0000_0002p+0, 0x1.8p+5, 0x4004_8000_0000_0000, 0x0000_0000_0000_0003);
108}
109
110const qnan80 = @bitCast(f80, @bitCast(u80, math.nan(f80)) | (1 << (math.floatFractionalBits(f80) - 1)));
111
112fn test__mulxf3(a: f80, b: f80, expected: u80) !void {
113 const x = __mulxf3(a, b);
114 const rep = @bitCast(u80, x);
115
116 if (rep == expected)
117 return;
118
119 if (math.isNan(@bitCast(f80, expected)) and math.isNan(x))
120 return; // We don't currently test NaN payload propagation
121
122 return error.TestFailed;
123}
124
125test "mulxf3" {
126 // NaN * any = NaN
127 try test__mulxf3(qnan80, 0x1.23456789abcdefp+5, @bitCast(u80, qnan80));
128 try test__mulxf3(@bitCast(f80, @as(u80, 0x7fff_8000_8000_3000_0000)), 0x1.23456789abcdefp+5, @bitCast(u80, qnan80));
129
130 // any * NaN = NaN
131 try test__mulxf3(0x1.23456789abcdefp+5, qnan80, @bitCast(u80, qnan80));
132 try test__mulxf3(0x1.23456789abcdefp+5, @bitCast(f80, @as(u80, 0x7fff_8000_8000_3000_0000)), @bitCast(u80, qnan80));
133
134 // NaN * inf = NaN
135 try test__mulxf3(qnan80, math.inf(f80), @bitCast(u80, qnan80));
136
137 // inf * NaN = NaN
138 try test__mulxf3(math.inf(f80), qnan80, @bitCast(u80, qnan80));
139
140 // inf * inf = inf
141 try test__mulxf3(math.inf(f80), math.inf(f80), @bitCast(u80, math.inf(f80)));
142
143 // inf * -inf = -inf
144 try test__mulxf3(math.inf(f80), -math.inf(f80), @bitCast(u80, -math.inf(f80)));
145
146 // -inf + inf = -inf
147 try test__mulxf3(-math.inf(f80), math.inf(f80), @bitCast(u80, -math.inf(f80)));
148
149 // inf * any = inf
150 try test__mulxf3(math.inf(f80), 0x1.2335653452436234723489432abcdefp+5, @bitCast(u80, math.inf(f80)));
151
152 // any * inf = inf
153 try test__mulxf3(0x1.2335653452436234723489432abcdefp+5, math.inf(f80), @bitCast(u80, math.inf(f80)));
154
155 // any * any
156 try test__mulxf3(0x1.0p+0, 0x1.dcba987654321p+5, 0x4004_ee5d_4c3b_2a19_0800);
157 try test__mulxf3(0x1.0000_0000_0000_0004p+0, 0x1.8p+5, 0x4004_C000_0000_0000_0003); // exact
158
159 try test__mulxf3(0x1.0000_0000_0000_0002p+0, 0x1.0p+5, 0x4004_8000_0000_0000_0001); // exact
160 try test__mulxf3(0x1.0000_0000_0000_0002p+0, 0x1.7ffep+5, 0x4004_BFFF_0000_0000_0001); // round down
161 try test__mulxf3(0x1.0000_0000_0000_0002p+0, 0x1.8p+5, 0x4004_C000_0000_0000_0002); // round up to even
162 try test__mulxf3(0x1.0000_0000_0000_0002p+0, 0x1.8002p+5, 0x4004_C001_0000_0000_0002); // round up
163 try test__mulxf3(0x1.0000_0000_0000_0002p+0, 0x1.0p+6, 0x4005_8000_0000_0000_0001); // exact
164
165 try test__mulxf3(0x1.0000_0001p+0, 0x1.0000_0001p+0, 0x3FFF_8000_0001_0000_0000); // round down to even
166 try test__mulxf3(0x1.0000_0001p+0, 0x1.0000_0001_0002p+0, 0x3FFF_8000_0001_0001_0001); // round up
100167}