Hopefully completely fix printf-family %a rounding (#1287)

The a conversion specifier to printf had some issues w.r.t. rounding, in
particular in edge cases w.r.t. "to nearest, ties to even" rounding (for
instance, "%.1a" with 0x1.78p+4 outputted 0x1.7p+4 instead of 0x1.8p+4).

This patch fixes this and adds several tests w.r.t ties to even rounding

Author: GabrielRavier

libc/stdio/fmt.c

@@ -565,12 +565,12 @@ static int __fmt_stoa(int out(const char *, void *, size_t), void *arg,
 static void __fmt_dfpbits(union U *u, struct FPBits *b) {
   int ex, i;
   b->fpi = kFpiDbl;
-  // Uncomment this if needed in the future - we currently do not need it, as
-  // the only reason we need it in __fmt_ldfpbits is because gdtoa reads
-  // fpi.rounding to determine rounding (which dtoa does not need as it directly
-  // reads FLT_ROUNDS)
-  // if (FLT_ROUNDS != -1)
-  //   b->fpi.rounding = FLT_ROUNDS;
+
+  // dtoa doesn't need this, unlike gdtoa, but we use it for __fmt_bround
+  i = FLT_ROUNDS;
+  if (i != -1)
+    b->fpi.rounding = i;
+
   b->sign = u->ui[1] & 0x80000000L;
   b->bits[1] = u->ui[1] & 0xfffff;
   b->bits[0] = u->ui[0];
@@ -616,10 +616,14 @@ static void __fmt_ldfpbits(union U *u, struct FPBits *b) {
 #error "unsupported architecture"
 #endif
   b->fpi = kFpiLdbl;
+
   // gdtoa doesn't check for FLT_ROUNDS but for fpi.rounding (which has the
   // same valid values as FLT_ROUNDS), so handle this here
-  if (FLT_ROUNDS != -1)
-    b->fpi.rounding = FLT_ROUNDS;
+  // (we also use this in __fmt_bround now)
+  i = FLT_ROUNDS;
+  if (i != -1)
+    b->fpi.rounding = i;
+
   b->sign = sex & 0x8000;
   if ((ex = sex & 0x7fff) != 0) {
     if (ex != 0x7fff) {
@@ -692,7 +696,6 @@ static int __fmt_bround(struct FPBits *b, int prec, int prec1) {
   uint32_t *bits, t;
   int i, j, k, m, n;
   bool inc = false;
-  int current_rounding_mode;
   m = prec1 - prec;
   bits = b->bits;
   k = m - 1;
@@ -701,22 +704,24 @@ static int __fmt_bround(struct FPBits *b, int prec, int prec1) {
   // always know in which direction we must round because of the current
   // rounding mode (note that if the correct value for inc is `false` then it
   // doesn't need to be set as we have already done so above)
-  // The last one handles rounding to nearest
-  current_rounding_mode = fegetround();
-  if (current_rounding_mode == FE_TOWARDZERO ||
-      (current_rounding_mode == FE_UPWARD && b->sign) ||
-      (current_rounding_mode == FE_DOWNWARD && !b->sign))
+  // They use the FLT_ROUNDS value, which are the same as gdtoa's FPI_Round_*
+  // enum values
+  if (b->fpi.rounding == FPI_Round_zero ||
+      (b->fpi.rounding == FPI_Round_up && b->sign) ||
+      (b->fpi.rounding == FPI_Round_down && !b->sign))
     goto have_inc;
-  if ((current_rounding_mode == FE_UPWARD && !b->sign) ||
-      (current_rounding_mode == FE_DOWNWARD && b->sign)) {
-    inc = true;
-    goto have_inc;
-  }
+  if ((b->fpi.rounding == FPI_Round_up && !b->sign) ||
+      (b->fpi.rounding == FPI_Round_down && b->sign))
+    goto inc_true;
 
   if ((t = bits[k >> 3] >> (j = (k & 7) * 4)) & 8) {
     if (t & 7)
       goto inc_true;
-    if (j && bits[k >> 3] << (32 - j))
+    // ((1 << (j * 4)) - 1) will mask appropriately for the lower bits
+    if ((bits[k >> 3] & ((1 << (j * 4)) - 1)) != 0)
+      goto inc_true;
+    // If exactly halfway and all lower bits are zero (tie), round to even
+    if ((bits[k >> 3] >> (j + 1) * 4) & 1)
       goto inc_true;
     while (k >= 8) {
       k -= 8;

test/libc/stdio/snprintf_test.c

@@ -217,39 +217,37 @@ TEST(snprintf, testLongDoubleRounding) {
   ASSERT_EQ(0, fesetround(previous_rounding));
 }
 
+void check_a_conversion_specifier_prec_1(const char *result_str, double value) {
+  char buf[30] = {0};
+  int i = snprintf(buf, sizeof(buf), "%.1a", value);
+
+  ASSERT_EQ(strlen(result_str), i);
+  ASSERT_STREQ(result_str, buf);
+}
+
 TEST(snprintf, testAConversionSpecifierRounding) {
   int previous_rounding = fegetround();
-  ASSERT_EQ(0, fesetround(FE_DOWNWARD));
 
-  char buf[20];
-  int i = snprintf(buf, sizeof(buf), "%.1a", 0x1.fffffp+4);
-  ASSERT_EQ(8, i);
-  ASSERT_STREQ("0x1.fp+4", buf);
+  ASSERT_EQ(0, fesetround(FE_DOWNWARD));
+  check_a_conversion_specifier_prec_1("0x1.fp+4", 0x1.fffffp+4);
 
   ASSERT_EQ(0, fesetround(FE_UPWARD));
-
-  i = snprintf(buf, sizeof(buf), "%.1a", 0x1.f8p+4);
-  ASSERT_EQ(8, i);
-  ASSERT_STREQ("0x2.0p+4", buf);
+  check_a_conversion_specifier_prec_1("0x2.0p+4", 0x1.f8p+4);
 
   ASSERT_EQ(0, fesetround(previous_rounding));
 }
 
-// This test currently fails because of rounding issues
-// If that ever gets fixed, uncomment this
-/*
+// This test specifically checks that we round to even, accordingly to IEEE
+// rules
 TEST(snprintf, testAConversionSpecifier) {
-  char buf[20];
-  int i = snprintf(buf, sizeof(buf), "%.1a", 0x1.7800000000001p+4);
-  ASSERT_EQ(8, i);
-  ASSERT_STREQ("0x1.8p+4", buf);
-
-  memset(buf, 0, sizeof(buf));
-  i = snprintf(buf, sizeof(buf), "%.1a", 0x1.78p+4);
-  ASSERT_EQ(8, i);
-  ASSERT_STREQ("0x1.8p+4", buf);
+  check_a_conversion_specifier_prec_1("0x1.8p+4", 0x1.7800000000001p+4);
+  check_a_conversion_specifier_prec_1("0x1.8p+4", 0x1.78p+4);
+  check_a_conversion_specifier_prec_1("0x1.8p+4", 0x1.88p+4);
+  check_a_conversion_specifier_prec_1("0x1.6p+4", 0x1.58p+4);
+  check_a_conversion_specifier_prec_1("0x1.6p+4", 0x1.68p+4);
+  check_a_conversion_specifier_prec_1("0x1.ap+4", 0x1.98p+4);
+  check_a_conversion_specifier_prec_1("0x1.ap+4", 0x1.a8p+4);
 }
-*/
 
 TEST(snprintf, apostropheFlag) {
   char buf[20];