authorgravatar for mrjbq7@gmail.comJohn Benediktsson <mrjbq7@gmail.com> 2026-10-12 13:50:27-07:00
committergravatar for alex@alexrp.comAlex Rønne Petersen <alex@alexrp.com> 2026-10-14 21:08:33+02:00
logcc958207734b68df15cd6ced350281753e97806c
treefa0d39edf693f3ff85b1aaf327e97218001f0027
parent4428d52ea145ed0f14af2e39569026c2a708e547
signature Signed by SSH key SHA256:7B/LJ7bpR1eX8aCXSr4mtd5M45VMPKcx9zY8e95b5QM

compiler_rt: fix fma rounding and add regression tests

1 files changed, 75 insertions(+), 6 deletions(-)

lib/compiler_rt/fma.zig+75-6
......@@ -8,6 +8,7 @@
88const std = @import("std");
99const math = std.math;
1010const expect = std.testing.expect;
11const expectEqual = std.testing.expectEqual;
1112const compiler_rt = @import("../compiler_rt.zig");
1213const symbol = compiler_rt.symbol;
1314
......@@ -40,8 +41,10 @@ pub fn fma_f32(x: f32, y: f32, z: f32) f32 {
4041 if ((u & 0x1FFFFFFF) != 0x10000000 or e == 0x7FF or (xy_z - xy == z and xy_z - z == xy)) {
4142 return @floatCast(xy_z);
4243 } else {
43 // TODO: Handle inexact case with double-rounding
44 return @floatCast(xy_z);
44 // Correct double rounding at an f32 midpoint.
45 const negative = u >> 63 != 0;
46 const err = if (negative == (z > xy)) xy - xy_z + z else z - xy_z + xy;
47 return @floatCast(@as(f64, @bitCast(if (negative == (err < 0)) u + 1 else u - 1)));
4548 }
4649}
4750
......@@ -74,6 +77,7 @@ pub fn fma_f64(x: f64, y: f64, z: f64) f64 {
7477 var zs = x3.significand;
7578
7679 var spread = ex + ey - ez;
80 if (spread < -53) return z;
7781 if (spread <= 53 * 2) {
7882 zs = math.scalbn(zs, -spread);
7983 } else {
......@@ -85,7 +89,7 @@ pub fn fma_f64(x: f64, y: f64, z: f64) f64 {
8589 spread = ex + ey;
8690
8791 if (r.hi == 0.0) {
88 return xy.hi + zs + math.scalbn(xy.lo, spread);
92 return if (xy.lo == 0.0) xy.hi + zs else math.scalbn(xy.lo, spread);
8993 }
9094
9195 const adj = add_adjusted(r.lo, xy.lo);
......@@ -139,6 +143,7 @@ pub fn fma_f128(x: f128, y: f128, z: f128) f128 {
139143 var zs = x3.significand;
140144
141145 var spread = ex + ey - ez;
146 if (spread < -113) return z;
142147 if (spread <= 113 * 2) {
143148 zs = math.scalbn(zs, -spread);
144149 } else {
......@@ -150,7 +155,7 @@ pub fn fma_f128(x: f128, y: f128, z: f128) f128 {
150155 spread = ex + ey;
151156
152157 if (r.hi == 0.0) {
153 return xy.hi + zs + math.scalbn(xy.lo, spread);
158 return if (xy.lo == 0.0) xy.hi + zs else math.scalbn(xy.lo, spread);
154159 }
155160
156161 const adj = add_adjusted128(r.lo, xy.lo);
......@@ -224,7 +229,7 @@ fn add_and_denorm(a: f64, b: f64, scale: i32) f64 {
224229 if (sum.lo != 0) {
225230 var uhii: u64 = @bitCast(sum.hi);
226231 const bits_lost = -@as(i32, @intCast((uhii >> 52) & 0x7FF)) - scale + 1;
227 if ((bits_lost != 1) == (uhii & 1 != 0)) {
232 if ((bits_lost != 1) != (uhii & 1 != 0)) {
228233 const uloi: u64 = @bitCast(sum.lo);
229234 uhii = uhii + 1 - (((uhii ^ uloi) >> 62) & 2);
230235 sum.hi = @bitCast(uhii);
......@@ -291,7 +296,7 @@ fn add_and_denorm128(a: f128, b: f128, scale: i32) f128 {
291296 if (sum.lo != 0) {
292297 var uhii: u128 = @bitCast(sum.hi);
293298 const bits_lost = -@as(i32, @intCast((uhii >> 112) & 0x7FFF)) - scale + 1;
294 if ((bits_lost != 1) == (uhii & 1 != 0)) {
299 if ((bits_lost != 1) != (uhii & 1 != 0)) {
295300 const uloi: u128 = @bitCast(sum.lo);
296301 uhii = uhii + 1 - (((uhii ^ uloi) >> 126) & 2);
297302 sum.hi = @bitCast(uhii);
......@@ -335,6 +340,9 @@ test "32" {
335340 try expect(math.approxEqAbs(f32, fma_f32(37.45, 5.0, 9.124), 196.374004, epsilon));
336341 try expect(math.approxEqAbs(f32, fma_f32(89.123, 5.0, 9.124), 454.739005, epsilon));
337342 try expect(math.approxEqAbs(f32, fma_f32(123123.234375, 5.0, 9.124), 615625.295875, epsilon));
343
344 // Avoid double rounding when the f64 sum lands on an f32 midpoint.
345 try expectEqual(0x3fc00001, @as(u32, @bitCast(fma_f32(0x1.000002p0, 1.5, -0x1p-60))));
338346}
339347
340348test "64" {
......@@ -347,6 +355,45 @@ test "64" {
347355 try expect(math.approxEqAbs(f64, fma_f64(37.45, 5.0, 9.124), 196.374, epsilon));
348356 try expect(math.approxEqAbs(f64, fma_f64(89.123, 5.0, 9.124), 454.739, epsilon));
349357 try expect(math.approxEqAbs(f64, fma_f64(123123.234375, 5.0, 9.124), 615625.295875, epsilon));
358
359 // Special operands.
360 try expect(math.isNan(fma_f64(math.inf(f64), 0.0, 1.0)));
361 try expect(math.isNan(fma_f64(math.inf(f64), 1.0, -math.inf(f64))));
362 try expectEqual(math.inf(f64), fma_f64(2.0, 3.0, math.inf(f64)));
363 try expectEqual(math.inf(f64), fma_f64(math.floatMax(f64), 2.0, 0.0));
364 try expect(math.isNan(fma_f64(math.nan(f64), 1.0, 1.0)));
365
366 // Addend dominates product.
367 try expectEqual(-1.438e308, fma_f64(-0.6, math.floatTrueMin(f64), -1.438e308));
368 try expectEqual(1.0, fma_f64(0x1p-600, 0x1p-600, 1.0));
369 try expectEqual(math.floatMax(f64), fma_f64(math.floatTrueMin(f64), math.floatTrueMin(f64), math.floatMax(f64)));
370
371 // Subnormal result is rounded once.
372 try expectEqual(0x000d2971831e2f43, @as(u64, @bitCast(fma_f64(
373 @bitCast(@as(u64, 0xbbbcae314558d732)),
374 @bitCast(@as(u64, 0x843d607d0e518063)),
375 @bitCast(@as(u64, 0x800000ba24d8f3ba)),
376 ))));
377 try expectEqual(0x000516a34ce113eb, @as(u64, @bitCast(fma_f64(
378 @bitCast(@as(u64, 0x2af4b8a2ad5294f0)),
379 @bitCast(@as(u64, 0x14ef6dd2b62857c2)),
380 @bitCast(@as(u64, 0x000000248ff787f4)),
381 ))));
382 try expectEqual(0x800bccc9d0e68a53, @as(u64, @bitCast(fma_f64(
383 @bitCast(@as(u64, 0x3e81368dcd4cb557)),
384 @bitCast(@as(u64, 0x817514ce659bdad8)),
385 @bitCast(@as(u64, 0x800075cba3b3043c)),
386 ))));
387
388 // Sign of zero.
389 try expectEqual(0x8000000000000000, @as(u64, @bitCast(fma_f64(
390 @bitCast(@as(u64, 0x2001526ce7a829e5)),
391 @bitCast(@as(u64, 0xa00c9ac305fd618a)),
392 @bitCast(@as(u64, 0x001ef7cb10a75631)),
393 ))));
394 try expect(!math.signbit(fma_f64(1.0, 1.0, -1.0)));
395 try expect(math.signbit(fma_f64(-math.floatTrueMin(f64), 0.5, 0.0)));
396 try expect(math.signbit(fma_f64(-0.0, 1.0, -0.0)));
350397}
351398
352399test "128" {
......@@ -359,4 +406,26 @@ test "128" {
359406 try expect(math.approxEqAbs(f128, fma_f128(37.45, 5.0, 9.124), 196.374, epsilon));
360407 try expect(math.approxEqAbs(f128, fma_f128(89.123, 5.0, 9.124), 454.739, epsilon));
361408 try expect(math.approxEqAbs(f128, fma_f128(123123.234375, 5.0, 9.124), 615625.295875, epsilon));
409
410 // Addend dominates product.
411 try expectEqual(-1.0e4932, fma_f128(-0.6, math.floatTrueMin(f128), -1.0e4932));
412 try expectEqual(0x8000f2eb9b403234570b698042da365d, @as(u128, @bitCast(fma_f128(
413 @bitCast(@as(u128, 0xa4d80b69452301dd5a3efab045ca0371)),
414 @bitCast(@as(u128, 0x13a5efea97a6ece0f0b24ffab5fcdcb4)),
415 @bitCast(@as(u128, 0x8000f2eb9b403234570b698042da365d)),
416 ))));
417
418 // Subnormal result is rounded once.
419 try expectEqual(0x0000f5ea6135ffda377d6f2322102a7b, @as(u128, @bitCast(fma_f128(
420 @bitCast(@as(u128, 0x2d83ef38da3dc9f883f14ecb887541ce)),
421 @bitCast(@as(u128, 0x123129a960f91ee507caa6a6702fc33e)),
422 @bitCast(@as(u128, 0x0000f5ea6135ffda377d6eff2506556b)),
423 ))));
424
425 // Sign of zero.
426 try expectEqual(0x80000000000000000000000000000000, @as(u128, @bitCast(fma_f128(
427 @bitCast(@as(u128, 0x1fd32f92b54d7cefc06df6f8521f1bc3)),
428 @bitCast(@as(u128, 0x202bb49b61ca93fdb382f78b6a2a32e2)),
429 @bitCast(@as(u128, 0x8000816f87a4a3a99ed5f1db8e4c0404)),
430 ))));
362431}