Fix printf-family long double prec/rounding issues (#1283)

GabrielRavier · web-flow · commit 4754f200ee50 · 2024-09-07T18:26:04.000-07:00
Currently, in cosmopolitan, there is no handling of the current rounding
mode for long double conversions, such that round-to-nearest gets always
used, regardless of the current rounding mode. %Le also improperly calls
gdtoa with a too small precision (which led to relatively similar bugs).

This patch fixes these issues, in particular by modifying the FPI object
passed to gdtoa such that it is modifiable (so that __fmt can adjust its
rounding field to correspond to FLT_ROUNDS (note that this is not needed
for dtoa, which checks FLT_ROUNDS directly)) and ors STRTOG_Neg into the
kind field in both of the __fmt_dfpbits and __fmt_ldfpbits functions, as
the gdtoa function also depends on it to be able to accurately round any
negative arguments. The change to kind also requires a few other changes
to make sure kind's upper bits (which include STRTOG_Neg) are masked off
when attempting to only examine the lower bits' value. Furthermore, this
patch also makes exactly one change in gdtoa, which appears to be needed
to fix rounding issues with FE_TOWARDZERO (this seems like a gdtoa bug).

The patch also adds a few tests for these issues, along with also taking
the opportunity to clean up some of the previous tests to do the asserts
in the right order (i.e. with the first argument as the expected result,
and the second one being used as the value that it is compared against).
diff --git a/libc/stdio/fmt.c b/libc/stdio/fmt.c
@@ -53,6 +53,7 @@
 #include "libc/math.h"
 #include "libc/mem/mem.h"
 #include "libc/mem/reverse.internal.h"
+#include "libc/runtime/fenv.h"
 #include "libc/runtime/internal.h"
 #include "libc/serialize.h"
 #include "libc/str/str.h"
@@ -90,7 +91,7 @@
 
 struct FPBits {
   uint32_t bits[4];
-  const FPI *fpi;
+  FPI fpi;
   int sign;
   int ex;  // exponent
   int kind;
@@ -563,7 +564,13 @@ 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;
+  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;
   b->sign = u->ui[1] & 0x80000000L;
   b->bits[1] = u->ui[1] & 0xfffff;
   b->bits[0] = u->ui[0];
@@ -582,6 +589,7 @@ static void __fmt_dfpbits(union U *u, struct FPBits *b) {
   } else {
     i = STRTOG_Zero;
   }
+  i |= signbit(u->d) ? STRTOG_Neg : 0;
   b->kind = i;
   b->ex = ex - (0x3ff + 52);
 }
@@ -607,7 +615,11 @@ static void __fmt_ldfpbits(union U *u, struct FPBits *b) {
 #else
 #error "unsupported architecture"
 #endif
-  b->fpi = &kFpiLdbl;
+  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;
   b->sign = sex & 0x8000;
   if ((ex = sex & 0x7fff) != 0) {
     if (ex != 0x7fff) {
@@ -626,6 +638,7 @@ static void __fmt_ldfpbits(union U *u, struct FPBits *b) {
   } else {
     i = STRTOG_Zero;
   }
+  i |= signbit(u->ld) ? STRTOG_Neg : 0;
   b->kind = i;
   b->ex = ex - (0x3fff + (LDBL_MANT_DIG - 1));
 #endif
@@ -636,9 +649,9 @@ static int __fmt_fpiprec(struct FPBits *b) {
   const FPI *fpi;
   int i, j, k, m;
   uint32_t *bits;
-  if (b->kind == STRTOG_Zero)
+  if ((b->kind & STRTOG_Retmask) == STRTOG_Zero)
     return (b->ex = 0);
-  fpi = b->fpi;
+  fpi = &b->fpi;
   bits = b->bits;
   for (k = (fpi->nbits - 1) >> 2; k > 0; --k) {
     if ((bits[k >> 3] >> 4 * (k & 7)) & 0xf) {
@@ -1163,8 +1176,8 @@ int __fmt(void *fn, void *arg, const char *format, va_list va, int *wrote) {
         } else {
           un.ld = va_arg(va, long double);
           __fmt_ldfpbits(&un, &fpb);
-          s = s0 =
-              gdtoa(fpb.fpi, fpb.ex, fpb.bits, &fpb.kind, 3, prec, &decpt, &se);
+          s = s0 = gdtoa(&fpb.fpi, fpb.ex, fpb.bits, &fpb.kind, 3, prec, &decpt,
+                         &se);
         }
         if (s0 == NULL)
           return -1;
@@ -1181,7 +1194,10 @@ int __fmt(void *fn, void *arg, const char *format, va_list va, int *wrote) {
           } else if (flags & FLAGS_SPACE) {
             *q++ = ' ';
           }
-          memcpy(q, kSpecialFloats[fpb.kind == STRTOG_NaN][d >= 'a'], 4);
+          memcpy(q,
+                 kSpecialFloats[(fpb.kind & STRTOG_Retmask) == STRTOG_NaN]
+                               [d >= 'a'],
+                 4);
           flags &= ~(FLAGS_PRECISION | FLAGS_PLUS | FLAGS_HASH | FLAGS_SPACE);
           prec = 0;
           rc = __fmt_stoa(out, arg, s, flags, prec, width, signbit, qchar);
@@ -1278,7 +1294,7 @@ int __fmt(void *fn, void *arg, const char *format, va_list va, int *wrote) {
         } else {
           un.ld = va_arg(va, long double);
           __fmt_ldfpbits(&un, &fpb);
-          s = s0 = gdtoa(fpb.fpi, fpb.ex, fpb.bits, &fpb.kind, prec ? 2 : 0,
+          s = s0 = gdtoa(&fpb.fpi, fpb.ex, fpb.bits, &fpb.kind, prec ? 2 : 0,
                          prec, &decpt, &se);
         }
         if (s0 == NULL)
@@ -1326,8 +1342,8 @@ int __fmt(void *fn, void *arg, const char *format, va_list va, int *wrote) {
         } else {
           un.ld = va_arg(va, long double);
           __fmt_ldfpbits(&un, &fpb);
-          s = s0 = gdtoa(fpb.fpi, fpb.ex, fpb.bits, &fpb.kind, prec ? 2 : 0,
-                         prec, &decpt, &se);
+          s = s0 = gdtoa(&fpb.fpi, fpb.ex, fpb.bits, &fpb.kind, prec ? 2 : 0,
+                         prec + 1, &decpt, &se);
         }
         if (s0 == NULL)
           return -1;
@@ -1411,7 +1427,8 @@ int __fmt(void *fn, void *arg, const char *format, va_list va, int *wrote) {
           un.d = va_arg(va, double);
           __fmt_dfpbits(&un, &fpb);
         }
-        if (fpb.kind == STRTOG_Infinite || fpb.kind == STRTOG_NaN) {
+        if ((fpb.kind & STRTOG_Retmask) == STRTOG_Infinite ||
+            (fpb.kind & STRTOG_Retmask) == STRTOG_NaN) {
           s0 = 0;
           goto FormatDecpt9999Or32768;
         }
@@ -1429,7 +1446,7 @@ int __fmt(void *fn, void *arg, const char *format, va_list va, int *wrote) {
             i1 /= 10;
           }
         }
-        if (fpb.sign /* && (sign || fpb.kind != STRTOG_Zero) */) {
+        if (fpb.sign /* && (sign || (fpb.kind & STRTOG_Retmask) != STRTOG_Zero) */) {
           sign = '-';
         }
         if ((width -= bw + 5) > 0) {
diff --git a/test/libc/stdio/snprintf_test.c b/test/libc/stdio/snprintf_test.c
@@ -16,6 +16,7 @@
 │ TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR             │
 │ PERFORMANCE OF THIS SOFTWARE.                                                │
 ╚─────────────────────────────────────────────────────────────────────────────*/
+#include "libc/runtime/fenv.h"
 #include "libc/stdio/stdio.h"
 #include "libc/str/str.h"
 #include "libc/testlib/testlib.h"
@@ -24,64 +25,64 @@ TEST(snprintf, testVeryLargePrecision) {
   char buf[512] = {};
   int i = snprintf(buf, sizeof(buf), "%.9999u", 10);
 
-  ASSERT_EQ(i, 9999);
-  ASSERT_EQ(strlen(buf), 511);
+  ASSERT_EQ(9999, i);
+  ASSERT_EQ(511, strlen(buf));
 }
 
 TEST(snprintf, testPlusFlagOnChar) {
   char buf[10] = {};
   int i = snprintf(buf, sizeof(buf), "%+c", '=');
 
-  ASSERT_EQ(i, 1);
-  ASSERT_STREQ(buf, "=");
+  ASSERT_EQ(1, i);
+  ASSERT_STREQ("=", buf);
 }
 
 TEST(snprintf, testInf) {
   char buf[10] = {};
   int i = snprintf(buf, sizeof(buf), "%f", 1.0 / 0.0);
 
-  ASSERT_EQ(i, 3);
-  ASSERT_STREQ(buf, "inf");
+  ASSERT_EQ(3, i);
+  ASSERT_STREQ("inf", buf);
 
   memset(buf, '\0', 4);
   i = snprintf(buf, sizeof(buf), "%Lf", 1.0L / 0.0L);
-  ASSERT_EQ(i, 3);
-  ASSERT_STREQ(buf, "inf");
+  ASSERT_EQ(3, i);
+  ASSERT_STREQ("inf", buf);
 
   memset(buf, '\0', 4);
   i = snprintf(buf, sizeof(buf), "%e", 1.0 / 0.0);
-  ASSERT_EQ(i, 3);
-  ASSERT_STREQ(buf, "inf");
+  ASSERT_EQ(3, i);
+  ASSERT_STREQ("inf", buf);
 
   memset(buf, '\0', 4);
   i = snprintf(buf, sizeof(buf), "%Le", 1.0L / 0.0L);
-  ASSERT_EQ(i, 3);
-  ASSERT_STREQ(buf, "inf");
+  ASSERT_EQ(3, i);
+  ASSERT_STREQ("inf", buf);
 
   memset(buf, '\0', 4);
   i = snprintf(buf, sizeof(buf), "%g", 1.0 / 0.0);
-  ASSERT_EQ(i, 3);
-  ASSERT_STREQ(buf, "inf");
+  ASSERT_EQ(3, i);
+  ASSERT_STREQ("inf", buf);
 
   memset(buf, '\0', 4);
   i = snprintf(buf, sizeof(buf), "%Lg", 1.0L / 0.0L);
-  ASSERT_EQ(i, 3);
-  ASSERT_STREQ(buf, "inf");
+  ASSERT_EQ(3, i);
+  ASSERT_STREQ("inf", buf);
 
   for (i = 4; i < 10; ++i)
-    ASSERT_EQ(buf[i], '\0');
+    ASSERT_EQ('\0', buf[i]);
 }
 
 TEST(snprintf, testUppercaseCConversionSpecifier) {
   char buf[10] = {};
   int i = snprintf(buf, sizeof(buf), "%C", L'a');
 
-  ASSERT_EQ(i, 1);
-  ASSERT_STREQ(buf, "a");
+  ASSERT_EQ(1, i);
+  ASSERT_STREQ("a", buf);
 
   i = snprintf(buf, sizeof(buf), "%C", L'☺');
-  ASSERT_EQ(i, 3);
-  ASSERT_STREQ(buf, "☺");
+  ASSERT_EQ(3, i);
+  ASSERT_STREQ("☺", buf);
 }
 
 // Make sure we don't va_arg the wrong argument size on wide character
@@ -92,96 +93,126 @@ TEST(snprintf,
   int i = snprintf(buf, sizeof(buf), "%d%d%d%d%d%d%d%d%lc%d", 0, 0, 0, 0, 0, 0,
                    0, 0, L'x', 1);
 
-  ASSERT_EQ(i, 10);
-  ASSERT_STREQ(buf, "00000000x1");
+  ASSERT_EQ(10, i);
+  ASSERT_STREQ("00000000x1", buf);
 
   memset(buf, 0, sizeof(buf));
   i = snprintf(buf, sizeof(buf), "%d%d%d%d%d%d%d%d%C%d", 0, 0, 0, 0, 0, 0, 0, 0,
                L'x', 1);
-  ASSERT_EQ(i, 10);
-  ASSERT_STREQ(buf, "00000000x1");
+  ASSERT_EQ(10, i);
+  ASSERT_STREQ("00000000x1", buf);
 }
 
 static void check_n_buffer_contents(char buf[350]) {
   for (int i = 0; i < 284; ++i)
-    ASSERT_EQ(buf[i], ' ');
-  ASSERT_STREQ(&buf[284], "428463");
+    ASSERT_EQ(' ', buf[i]);
+  ASSERT_STREQ("428463", &buf[284]);
   for (int i = 290; i < 350; ++i)
-    ASSERT_EQ(buf[i], '\0');
+    ASSERT_EQ('\0', buf[i]);
 }
 
 TEST(snprintf, testNConversionSpecifier) {
   char buf[350] = {};
 
   int n_res_int = -1;
   int i = snprintf(buf, sizeof(buf), "%286d%d%n%d", 42, 84, &n_res_int, 63);
-  ASSERT_EQ(i, 290);
+  ASSERT_EQ(290, i);
   check_n_buffer_contents(buf);
-  ASSERT_EQ(n_res_int, 288);
+  ASSERT_EQ(288, n_res_int);
 
   memset(&buf, '\0', sizeof(buf));
   long n_res_long = -1;
   i = snprintf(buf, sizeof(buf), "%286ld%ld%ln%ld", 42L, 84L, &n_res_long, 63L);
-  ASSERT_EQ(i, 290);
+  ASSERT_EQ(290, i);
   check_n_buffer_contents(buf);
-  ASSERT_EQ(n_res_long, 288);
+  ASSERT_EQ(288, n_res_long);
 
   memset(&buf, '\0', sizeof(buf));
   long long n_res_long_long = -1;
   i = snprintf(buf, sizeof(buf), "%286lld%lld%lln%lld", 42LL, 84LL,
                &n_res_long_long, 63LL);
-  ASSERT_EQ(i, 290);
+  ASSERT_EQ(290, i);
   check_n_buffer_contents(buf);
-  ASSERT_EQ(n_res_long_long, 288);
+  ASSERT_EQ(288, n_res_long_long);
 
   ASSERT_EQ(sizeof(short), 2);
   ASSERT_EQ(sizeof(int), 4);
   memset(&buf, '\0', sizeof(buf));
   short n_res_short = -1;
   i = snprintf(buf, sizeof(buf), "%286hd%hd%hn%hd", (42 | 0xFFFF0000),
                (84 | 0xFFFF0000), &n_res_short, (63 | 0xFFFF0000));
-  ASSERT_EQ(i, 290);
+  ASSERT_EQ(290, i);
   check_n_buffer_contents(buf);
-  ASSERT_EQ(n_res_short, 288);
+  ASSERT_EQ(288, n_res_short);
 
   ASSERT_EQ(sizeof(unsigned char), 1);
   memset(&buf, '\0', sizeof(buf));
   signed char n_res_char = -1;
   i = snprintf(buf, sizeof(buf), "%286hhd%hhd%hhn%hhd", (42 | 0xFFFFFF00),
                (84 | 0xFFFFFF00), &n_res_char, (63 | 0xFFFFFF00));
-  ASSERT_EQ(i, 290);
+  ASSERT_EQ(290, i);
   check_n_buffer_contents(buf);
-  ASSERT_EQ(n_res_char, (signed char)288);
+  ASSERT_EQ((signed char)288, n_res_char);
 
   memset(&buf, '\0', sizeof(buf));
   ssize_t n_res_size_t = -1;
   i = snprintf(buf, sizeof(buf), "%286zd%zd%zn%zd", (ssize_t)42, (ssize_t)84,
                &n_res_size_t, (ssize_t)63);
-  ASSERT_EQ(i, 290);
+  ASSERT_EQ(290, i);
   check_n_buffer_contents(buf);
-  ASSERT_EQ(n_res_size_t, 288);
+  ASSERT_EQ(288, n_res_size_t);
 
   memset(&buf, '\0', sizeof(buf));
   intmax_t n_res_intmax_t = -1;
   i = snprintf(buf, sizeof(buf), "%286jd%jd%jn%jd", (intmax_t)42, (intmax_t)84,
                &n_res_intmax_t, (intmax_t)63);
-  ASSERT_EQ(i, 290);
+  ASSERT_EQ(290, i);
   check_n_buffer_contents(buf);
-  ASSERT_EQ(n_res_intmax_t, 288);
+  ASSERT_EQ(288, n_res_intmax_t);
 
   memset(&buf, '\0', sizeof(buf));
   int128_t n_res_int128_t = -1;
   i = snprintf(buf, sizeof(buf), "%286jjd%jjd%jjn%jjd", (int128_t)42,
                (int128_t)84, &n_res_int128_t, (int128_t)63);
-  ASSERT_EQ(i, 290);
+  ASSERT_EQ(290, i);
   check_n_buffer_contents(buf);
-  ASSERT_EQ(n_res_int128_t, 288);
+  ASSERT_EQ(288, n_res_int128_t);
 
   memset(&buf, '\0', sizeof(buf));
   ptrdiff_t n_res_ptrdiff_t = -1;
   i = snprintf(buf, sizeof(buf), "%286td%td%tn%td", (ptrdiff_t)42,
                (ptrdiff_t)84, &n_res_ptrdiff_t, (ptrdiff_t)63);
-  ASSERT_EQ(i, 290);
+  ASSERT_EQ(290, i);
   check_n_buffer_contents(buf);
-  ASSERT_EQ(n_res_ptrdiff_t, 288);
+  ASSERT_EQ(288, n_res_ptrdiff_t);
+}
+
+TEST(snprintf, testLongDoubleEConversionSpecifier) {
+  char buf[20] = {};
+  int i = snprintf(buf, sizeof(buf), "%Le", 1234567.8L);
+
+  ASSERT_EQ(12, i);
+  ASSERT_STREQ("1.234568e+06", buf);
+}
+
+TEST(snprintf, testLongDoubleRounding) {
+  int previous_rounding = fegetround();
+  ASSERT_EQ(0, fesetround(FE_DOWNWARD));
+
+  char buf[20];
+  int i = snprintf(buf, sizeof(buf), "%.3Lf", 4.4375L);
+  ASSERT_EQ(5, i);
+  ASSERT_STREQ("4.437", buf);
+
+  i = snprintf(buf, sizeof(buf), "%.3Lf", -4.4375L);
+  ASSERT_EQ(6, i);
+  ASSERT_STREQ("-4.438", buf);
+
+  ASSERT_EQ(0, fesetround(FE_TOWARDZERO));
+
+  i = snprintf(buf, sizeof(buf), "%.3Lf", -4.4375L);
+  ASSERT_EQ(6, i);
+  ASSERT_STREQ("-4.437", buf);
+
+  ASSERT_EQ(0, fesetround(previous_rounding));
 }
diff --git a/third_party/gdtoa/gdtoa.c b/third_party/gdtoa/gdtoa.c
@@ -293,7 +293,9 @@ gdtoa(const FPI *fpi, int be, ULong *bits, int *kindp, int mode, int ndigits, in
 	else if ( (rdir = fpi->rounding - 1) !=0) {
 		if (rdir < 0)
 			rdir = 2;
-		if (kind & STRTOG_Neg)
+		// note that we check for fpi->rounding == 0 as in that case we
+		// must *always* round towards 0, i.e. downwards, with rdir = 2
+		if (kind & STRTOG_Neg && fpi->rounding != 0)
 			rdir = 3 - rdir;
 	}
 	/* Now rdir = 0 ==> round near, 1 ==> round up, 2 ==> round down. */
