Merge pull request #6722 from dhalbert/micropython-float-print-fix

tannewt · web-flow · commit 741a5c2beca8 · 2022-08-10T09:32:28.000-07:00
py/formatfloat: Format all whole-number floats exactly.
diff --git a/py/formatfloat.c b/py/formatfloat.c
@@ -25,6 +25,7 @@
  */
 
 #include "py/mpconfig.h"
+#include "py/misc.h"
 #if MICROPY_FLOAT_IMPL != MICROPY_FLOAT_IMPL_NONE
 
 #include <assert.h>
@@ -96,7 +97,16 @@ static inline int fp_isless1(float x) {
 #define fp_iszero(x) (x == 0)
 #define fp_isless1(x) (x < 1.0)
 
-#endif
+#endif // MICROPY_FLOAT_IMPL == MICROPY_FLOAT_IMPL_FLOAT/DOUBLE
+
+static inline int fp_ge_eps(FPTYPE x, FPTYPE y) {
+    mp_float_union_t fb_y = {y};
+    // Back off 2 eps.
+    // This is valid for almost all values, but in practice
+    // it's only used when y = 1eX for X>=0.
+    fb_y.i -= 2;
+    return x >= fb_y.f;
+}
 
 static const FPTYPE g_pos_pow[] = {
     #if FPDECEXP > 32
@@ -173,6 +183,7 @@ int mp_format_float(FPTYPE f, char *buf, size_t buf_size, char fmt, int prec, ch
     int num_digits = 0;
     const FPTYPE *pos_pow = g_pos_pow;
     const FPTYPE *neg_pow = g_neg_pow;
+    int signed_e = 0;
 
     if (fp_iszero(f)) {
         e = 0;
@@ -192,40 +203,32 @@ int mp_format_float(FPTYPE f, char *buf, size_t buf_size, char fmt, int prec, ch
             }
         }
     } else if (fp_isless1(f)) {
-        // We need to figure out what an integer digit will be used
-        // in case 'f' is used (or we revert other format to it below).
-        // As we just tested number to be <1, this is obviously 0,
-        // but we can round it up to 1 below.
-        char first_dig = '0';
-        if (f >= FPROUND_TO_ONE) {
-            first_dig = '1';
-        }
-
+        FPTYPE f_mod = f;
         // Build negative exponent
         for (e = 0, e1 = FPDECEXP; e1; e1 >>= 1, pos_pow++, neg_pow++) {
-            if (*neg_pow > f) {
+            if (*neg_pow > f_mod) {
                 e += e1;
-                f *= *pos_pow;
+                f_mod *= *pos_pow;
             }
         }
+
         char e_sign_char = '-';
-        if (fp_isless1(f) && f >= FPROUND_TO_ONE) {
-            f = FPCONST(1.0);
+        if (fp_isless1(f_mod) && f_mod >= FPROUND_TO_ONE) {
+            f_mod = FPCONST(1.0);
             if (e == 0) {
                 e_sign_char = '+';
             }
-        } else if (fp_isless1(f)) {
+        } else if (fp_isless1(f_mod)) {
             e++;
-            f *= FPCONST(10.0);
+            f_mod *= FPCONST(10.0);
         }
 
         // If the user specified 'g' format, and e is <= 4, then we'll switch
         // to the fixed format ('f')
 
         if (fmt == 'f' || (fmt == 'g' && e <= 4)) {
             fmt = 'f';
-            dec = -1;
-            *s++ = first_dig;
+            dec = 0;
 
             if (org_fmt == 'g') {
                 prec += (e - 1);
@@ -237,13 +240,8 @@ int mp_format_float(FPTYPE f, char *buf, size_t buf_size, char fmt, int prec, ch
             }
 
             num_digits = prec;
-            if (num_digits) {
-                *s++ = '.';
-                while (--e && num_digits) {
-                    *s++ = '0';
-                    num_digits--;
-                }
-            }
+            signed_e = 0;
+            ++num_digits;
         } else {
             // For e & g formats, we'll be printing the exponent, so set the
             // sign.
@@ -256,22 +254,29 @@ int mp_format_float(FPTYPE f, char *buf, size_t buf_size, char fmt, int prec, ch
                     prec++;
                 }
             }
+            signed_e = -e;
         }
     } else {
-        // Build positive exponent
-        for (e = 0, e1 = FPDECEXP; e1; e1 >>= 1, pos_pow++, neg_pow++) {
-            if (*pos_pow <= f) {
+        // Build positive exponent.
+        // We don't modify f at this point to avoid innaccuracies from
+        // scaling it.  Instead, we find the product of powers of 10
+        // that is not greater than it, and use that to start the
+        // mantissa.
+        FPTYPE u_base = FPCONST(1.0);
+        for (e = 0, e1 = FPDECEXP; e1; e1 >>= 1, pos_pow++) {
+            FPTYPE next_u = u_base * *pos_pow;
+            // fp_ge_eps performs "f >= (next_u - 2eps)" so that if, for
+            // numerical reasons, f is very close to a power of ten but
+            // not strictly equal, we still treat it as that power of 10.
+            // The comparison was failing for maybe 10% of 1eX values, but
+            // although rounding fixed many of them, there were still some
+            // rendering as 9.99999998e(X-1).
+            if (fp_ge_eps(f, next_u)) {
+                u_base = next_u;
                 e += e1;
-                f *= *neg_pow;
             }
         }
 
-        // It can be that f was right on the edge of an entry in pos_pow needs to be reduced
-        if ((int)f >= 10) {
-            e += 1;
-            f *= FPCONST(0.1);
-        }
-
         // If the user specified fixed format (fmt == 'f') and e makes the
         // number too big to fit into the available buffer, then we'll
         // switch to the 'e' format.
@@ -310,15 +315,15 @@ int mp_format_float(FPTYPE f, char *buf, size_t buf_size, char fmt, int prec, ch
         } else {
             e_sign = '+';
         }
+        signed_e = e;
     }
     if (prec < 0) {
         // This can happen when the prec is trimmed to prevent buffer overflow
         prec = 0;
     }
 
-    // We now have num.f as a floating point number between >= 1 and < 10
-    // (or equal to zero), and e contains the absolute value of the power of
-    // 10 exponent. and (dec + 1) == the number of dgits before the decimal.
+    // At this point e contains the absolute value of the power of 10 exponent.
+    // (dec + 1) == the number of dgits before the decimal.
 
     // For e, prec is # digits after the decimal
     // For f, prec is # digits after the decimal
@@ -336,25 +341,63 @@ int mp_format_float(FPTYPE f, char *buf, size_t buf_size, char fmt, int prec, ch
         num_digits = prec;
     }
 
-    // Print the digits of the mantissa
-    for (int i = 0; i < num_digits; ++i, --dec) {
-        int32_t d = (int32_t)f;
-        if (d < 0) {
-            *s++ = '0';
-        } else {
+    if (signed_e < 0) {
+        // The algorithm below treats numbers smaller than 1 by scaling them
+        // repeatedly by 10 to bring the new digit to the top.  Our input number
+        // was smaller than 1, so scale it up to be 1 <= f < 10.
+        FPTYPE u_base = FPCONST(1.0);
+        const FPTYPE *pow_u = g_pos_pow;
+        for (int m = FPDECEXP; m; m >>= 1, pow_u++) {
+            if (m & e) {
+                u_base *= *pow_u;
+            }
+        }
+        f *= u_base;
+    }
+
+    int d = 0;
+    int num_digits_left = num_digits;
+    for (int digit_index = signed_e; num_digits_left >= 0; --digit_index) {
+        FPTYPE u_base = FPCONST(1.0);
+        if (digit_index > 0) {
+            // Generate 10^digit_index for positive digit_index.
+            const FPTYPE *pow_u = g_pos_pow;
+            int target_index = digit_index;
+            for (int m = FPDECEXP; m; m >>= 1, pow_u++) {
+                if (m & target_index) {
+                    u_base *= *pow_u;
+                }
+            }
+        }
+        for (d = 0; d < 9; ++d) {
+            // This is essentially "if (f < u_base)", but with 2eps margin
+            // so that if f is just a tiny bit smaller, we treat it as
+            // equal (and accept the additional digit value).
+            if (!fp_ge_eps(f, u_base)) {
+                break;
+            }
+            f -= u_base;
+        }
+        // We calculate one more digit than we display, to use in rounding
+        // below.  So only emit the digit if it's one that we display.
+        if (num_digits_left > 0) {
+            // Emit this number (the leading digit).
             *s++ = '0' + d;
+            if (dec == 0 && prec > 0) {
+                *s++ = '.';
+            }
         }
-        if (dec == 0 && prec > 0) {
-            *s++ = '.';
+        --dec;
+        --num_digits_left;
+        if (digit_index <= 0) {
+            // Once we get below 1.0, we scale up f instead of calculting
+            // negative powers of 10 in u_base.  This provides better
+            // renditions of exact decimals like 1/16 etc.
+            f *= FPCONST(10.0);
         }
-        f -= (FPTYPE)d;
-        f *= FPCONST(10.0);
     }
-
-    // Round
-    // If we print non-exponential format (i.e. 'f'), but a digit we're going
-    // to round by (e) is too far away, then there's nothing to round.
-    if ((org_fmt != 'f' || e <= num_digits) && f >= FPCONST(5.0)) {
+    // Rounding.  If the next digit to print is >= 5, round up.
+    if (d >= 5) {
         char *rs = s;
         rs--;
         while (1) {
@@ -394,7 +437,10 @@ int mp_format_float(FPTYPE f, char *buf, size_t buf_size, char fmt, int prec, ch
                 }
             } else {
                 // Need at extra digit at the end to make room for the leading '1'
-                s++;
+                // but if we're at the buffer size limit, just drop the final digit.
+                if ((size_t)(s + 1 - buf) < buf_size) {
+                    s++;
+                }
             }
             char *ss = s;
             while (ss > rs) {
diff --git a/tests/float/float_format_ftoe.py b/tests/float/float_format_ftoe.py
@@ -0,0 +1,4 @@
+# check a case where rounding was suppressed inappropriately when "f" was
+# promoted to "e" for large numbers.
+v = 8.888e32
+print("%.2f" % v)  # '%.2f' format with e32 becomes '%.2e', expect 8.89e+32.
diff --git a/tests/float/float_format_ftoe.py.exp b/tests/float/float_format_ftoe.py.exp
@@ -0,0 +1 @@
+8.89e+32
diff --git a/tests/float/float_format_ints.py b/tests/float/float_format_ints.py
@@ -0,0 +1,31 @@
+# Test that integers format to exact values.
+
+for b in [13, 123, 457, 23456]:
+    for r in range(1, 10):
+        e_fmt = "{:." + str(r) + "e}"
+        f_fmt = "{:." + str(r) + "f}"
+        g_fmt = "{:." + str(r) + "g}"
+        for e in range(0, 5):
+            f = b * (10**e)
+            title = str(b) + " x 10^" + str(e)
+            print(title, "with format", e_fmt, "gives", e_fmt.format(f))
+            print(title, "with format", f_fmt, "gives", f_fmt.format(f))
+            print(title, "with format", g_fmt, "gives", g_fmt.format(f))
+
+# Check that powers of 10 (that fit in float32) format correctly.
+for i in range(31):
+    # It works to 12 digits on all platforms *except* qemu-arm, where
+    # 10^11 comes out as 10000000820 or something.
+    print("{:.7g}".format(float("1e" + str(i))))
+
+# 16777215 is 2^24 - 1, the largest integer that can be completely held
+# in a float32.
+print("{:f}".format(16777215))
+# 4294967040 = 16777215 * 128 is the largest integer that is exactly
+# represented by a float32 and that will also fit within a (signed) int32.
+# The upper bound of our integer-handling code is actually double this,
+# but that constant might cause trouble on systems using 32 bit ints.
+print("{:f}".format(2147483520))
+# Very large positive integers can be a test for precision and resolution.
+# This is a weird way to represent 1e38 (largest power of 10 for float32).
+print("{:.6e}".format(float("9" * 30 + "e8")))
diff --git a/tests/float/float_format_ints_doubleprec.py b/tests/float/float_format_ints_doubleprec.py
@@ -0,0 +1,15 @@
+# Test formatting of very large ints.
+# Relies on double-precision floats.
+
+import array
+import sys
+
+# Challenging way to express 1e200 and 1e100.
+print("{:.12e}".format(float("9" * 400 + "e-200")))
+print("{:.12e}".format(float("9" * 400 + "e-300")))
+
+# These correspond to the binary representation of 1e200 in float64s:
+v1 = 0x54B249AD2594C37D  # 1e100
+v2 = 0x6974E718D7D7625A  # 1e200
+print("{:.12e}".format(array.array("d", v1.to_bytes(8, sys.byteorder))[0]))
+print("{:.12e}".format(array.array("d", v2.to_bytes(8, sys.byteorder))[0]))
diff --git a/tests/run-tests.py b/tests/run-tests.py
@@ -426,6 +426,7 @@ def run_tests(pyb, tests, args, result_dir, num_threads=1):
     if upy_float_precision < 64:
         skip_tests.add("float/float_divmod.py")  # tested by float/float_divmod_relaxed.py instead
         skip_tests.add("float/float2int_doubleprec_intbig.py")
+        skip_tests.add("float/float_format_ints_doubleprec.py")
         skip_tests.add("float/float_parse_doubleprec.py")
 
     if not has_complex:
diff --git a/tools/tinytest-codegen.py b/tools/tinytest-codegen.py
@@ -100,6 +100,7 @@ def script_to_map(test_file):
     "float/float_divmod.py",
     # requires double precision floating point to work
     "float/float2int_doubleprec_intbig.py",
+    "float/float_format_ints_doubleprec.py",
     "float/float_parse_doubleprec.py",
     # inline asm FP tests (require Cortex-M4)
     "inlineasm/asmfpaddsub.py",
