authorgravatar for topolarity@tapscott.meCody Tapscott <topolarity@tapscott.me> 2022-04-18 20:23:02-07:00
committergravatar for topolarity@tapscott.meCody Tapscott <topolarity@tapscott.me> 2022-04-18 21:45:46-07:00
logb2950866b18916fe5e6e458b62ec294244fa2c2f
treeb526acac45878f23ce71f8351e5ccc02ab69ee38
parentd760cae2b1f77c239cc4b47a54124ff761aa8a7a

compiler_rt: Fix rounding/NaN handling for f80 add/sub

There were a few minor bugs in the rounding behavior and Inf/NaN handling for the f80 __addxf3 and __subtf3 functions. This change updates the original generic implementation to correctly handle f80 floats, including the explicit integer bit.

2 files changed, 107 insertions(+), 192 deletions(-)

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}