Merge PR-1684 from sverker/float_to_list-rounding-bug/OTP-14890

Fix rounding bug in float_to_list/2

2 files changed, 118 insertions, 57 deletions

diff --git a/erts/emulator/sys/common/erl_sys_common_misc.c b/erts/emulator/sys/common/erl_sys_common_misc.c
index 79f87eb3a9..ba03a11405 100644
--- a/erts/emulator/sys/common/erl_sys_common_misc.c
+++ b/erts/emulator/sys/common/erl_sys_common_misc.c
@@ -158,59 +158,50 @@ int
 sys_double_to_chars_fast(double f, char *buffer, int buffer_size, int decimals,
 			 int compact)
 {
-    /* Note that some C compilers don't support "static const" propagation
-     * so we use a defines */
-    #define SYS_DOUBLE_RND_CONST 0.55555555555555555
+    #define SYS_DOUBLE_RND_CONST 0.5
     #define FRAC_SIZE            52
     #define EXP_SIZE             11
-    #define EXP_MASK             ((1ll << EXP_SIZE) - 1)
+    #define EXP_MASK             (((Uint64)1 << EXP_SIZE) - 1)
     #define MAX_DECIMALS         (sizeof(cs_sys_double_pow10) \
 				   / sizeof(cs_sys_double_pow10[0]))
-    #define FRAC_MASK            ((1ll << FRAC_SIZE) - 1)
-    #define FRAC_MASK2           ((1ll << (FRAC_SIZE + 1)) - 1)
-    #define MAX_FLOAT            (1ll << (FRAC_SIZE+1))
+    #define FRAC_MASK            (((Uint64)1 << FRAC_SIZE) - 1)
+    #define FRAC_MASK2           (((Uint64)1 << (FRAC_SIZE + 1)) - 1)
+    #define MAX_FLOAT            ((Uint64)1 << (FRAC_SIZE+1))
 
     static const double cs_sys_double_pow10[] = {
-        SYS_DOUBLE_RND_CONST / 1ll,
-        SYS_DOUBLE_RND_CONST / 10ll,
-        SYS_DOUBLE_RND_CONST / 100ll,
-        SYS_DOUBLE_RND_CONST / 1000ll,
-        SYS_DOUBLE_RND_CONST / 10000ll,
-        SYS_DOUBLE_RND_CONST / 100000ll,
-        SYS_DOUBLE_RND_CONST / 1000000ll,
-        SYS_DOUBLE_RND_CONST / 10000000ll,
-        SYS_DOUBLE_RND_CONST / 100000000ll,
-        SYS_DOUBLE_RND_CONST / 1000000000ll,
-        SYS_DOUBLE_RND_CONST / 10000000000ll,
-        SYS_DOUBLE_RND_CONST / 100000000000ll,
-        SYS_DOUBLE_RND_CONST / 1000000000000ll,
-        SYS_DOUBLE_RND_CONST / 10000000000000ll,
-        SYS_DOUBLE_RND_CONST / 100000000000000ll,
-        SYS_DOUBLE_RND_CONST / 1000000000000000ll,
-        SYS_DOUBLE_RND_CONST / 10000000000000000ll,
-        SYS_DOUBLE_RND_CONST / 100000000000000000ll,
-        SYS_DOUBLE_RND_CONST / 1000000000000000000ll
+        SYS_DOUBLE_RND_CONST / 1e0,
+        SYS_DOUBLE_RND_CONST / 1e1,
+        SYS_DOUBLE_RND_CONST / 1e2,
+        SYS_DOUBLE_RND_CONST / 1e3,
+        SYS_DOUBLE_RND_CONST / 1e4,
+        SYS_DOUBLE_RND_CONST / 1e5,
+        SYS_DOUBLE_RND_CONST / 1e6,
+        SYS_DOUBLE_RND_CONST / 1e7,
+        SYS_DOUBLE_RND_CONST / 1e8,
+        SYS_DOUBLE_RND_CONST / 1e9,
+        SYS_DOUBLE_RND_CONST / 1e10,
+        SYS_DOUBLE_RND_CONST / 1e11,
+        SYS_DOUBLE_RND_CONST / 1e12,
+        SYS_DOUBLE_RND_CONST / 1e13,
+        SYS_DOUBLE_RND_CONST / 1e14,
+        SYS_DOUBLE_RND_CONST / 1e15,
+        SYS_DOUBLE_RND_CONST / 1e16,
+        SYS_DOUBLE_RND_CONST / 1e17,
+        SYS_DOUBLE_RND_CONST / 1e18
     };
 
-    long long mantissa, int_part = 0, frac_part = 0;
-    short exp;
+    Uint64 mantissa, int_part, frac_part;
+    int exp;
+    int fbits;
     int max;
     int neg;
     double fr;
-    union { long long L; double F; } x;
+    union { Uint64 L; double F; } x;
     char *p = buffer;
 
     if (decimals < 0)
         return -1;
 
-    /* Round the number to given decimal places. The number of 5's in the
-     * SYS_DOUBLE_RND_CONST constant is chosen such that adding any more 5's doesn't
-     * change the double precision of the number, i.e.:
-     * 1> term_to_binary(0.55555555555555555, [{minor_version, 1}]).
-     * <<131,70,63,225,199,28,113,199,28,114>>
-     * 2> term_to_binary(0.5555555555555555555, [{minor_version, 1}]).
-     * <<131,70,63,225,199,28,113,199,28,114>>
-     */
     if (f >= 0) {
         neg = 0;
         fr  = decimals < MAX_DECIMALS ? (f + cs_sys_double_pow10[decimals]) : f;
@@ -241,7 +232,7 @@ sys_double_to_chars_fast(double f, char *buffer, int buffer_size, int decimals,
     }
 
     exp      -= EXP_MASK >> 1;
-    mantissa |= (1ll << FRAC_SIZE);
+    mantissa |= ((Uint64)1 << FRAC_SIZE);
 
     /* Don't bother with optimizing too large numbers or too large precision */
     if (x.F > MAX_FLOAT || decimals >= MAX_DECIMALS) {
@@ -256,11 +247,16 @@ sys_double_to_chars_fast(double f, char *buffer, int buffer_size, int decimals,
         return p - buffer;
     } else if (exp >= FRAC_SIZE) {
         int_part  = mantissa << (exp - FRAC_SIZE);
+        frac_part = 0;
+        fbits = FRAC_SIZE;  /* not important as frac_part==0 */
     } else if (exp >= 0) {
-        int_part  = mantissa >> (FRAC_SIZE - exp);
-        frac_part = (mantissa << (exp + 1)) & FRAC_MASK2;
+        fbits = FRAC_SIZE - exp;
+        int_part  = mantissa >> fbits;
+        frac_part = mantissa & (((Uint64)1 << fbits) -1);
     } else /* if (exp < 0) */ {
-        frac_part = (mantissa & FRAC_MASK2) >> -(exp + 1);
+        int_part = 0;
+        frac_part = mantissa;
+        fbits = FRAC_SIZE - exp;
     }
 
     if (!int_part) {
@@ -270,9 +266,8 @@ sys_double_to_chars_fast(double f, char *buffer, int buffer_size, int decimals,
     } else {
         int ret, i, n;
         while (int_part != 0) {
-            long long j = int_part / 10;
-            *p++ = (char)(int_part - ((j << 3) + (j << 1)) + '0');
-            int_part = j;
+            *p++ = (char)((int_part % 10) + '0');
+            int_part /= 10;
         }
         if (neg)
             *p++ = '-';
@@ -298,11 +293,22 @@ sys_double_to_chars_fast(double f, char *buffer, int buffer_size, int decimals,
         max = decimals;
 
         for (i = 0; i < max; i++) {
-            /* frac_part *= 10; */
-            frac_part = (frac_part << 3) + (frac_part << 1);
-
-            *p++ = (char)((frac_part >> (FRAC_SIZE + 1)) + '0');
-            frac_part &= FRAC_MASK2;
+            if (frac_part > (ERTS_UINT64_MAX/5)) {
+                frac_part >>= 3;
+                fbits -= 3;
+            }
+
+            /* Multiply by 10 (5*2) to extract decimal digit as integer part */
+            frac_part *= 5;
+            fbits--;
+
+            if (fbits >= 64) {
+                *p++ = '0';
+            }
+            else {
+                *p++ = (char)((frac_part >> fbits) + '0');
+                frac_part &= ((Uint64)1 << fbits) - 1;
+            }
         }
 
         /* Delete trailing zeroes */
diff --git a/erts/emulator/test/num_bif_SUITE.erl b/erts/emulator/test/num_bif_SUITE.erl
index ac1a667185..d950179882 100644
--- a/erts/emulator/test/num_bif_SUITE.erl
+++ b/erts/emulator/test/num_bif_SUITE.erl
@@ -118,6 +118,7 @@ t_float(Config) when is_list(Config) ->
 %% Tests float_to_list/1, float_to_list/2, float_to_binary/1, float_to_binary/2
 
 t_float_to_string(Config) when is_list(Config) ->
+    rand_seed(),
     test_fts("0.00000000000000000000e+00", 0.0),
     test_fts("2.50000000000000000000e+01", 25.0),
     test_fts("2.50000000000000000000e+00", 2.5),
@@ -167,8 +168,8 @@ t_float_to_string(Config) when is_list(Config) ->
     test_fts("1.12300",1.123, [{decimals, 5}]),
     test_fts("1.123",1.123, [{decimals, 5}, compact]),
     test_fts("1.1234",1.1234,[{decimals, 6}, compact]),
-    test_fts("1.01",1.005, [{decimals, 2}]),
-    test_fts("-1.01",-1.005,[{decimals, 2}]),
+    test_fts("1.00",1.005, [{decimals, 2}]),  %% 1.005 is really 1.0049999999...
+    test_fts("-1.00",-1.005,[{decimals, 2}]),
     test_fts("0.999",0.999, [{decimals, 3}]),
     test_fts("-0.999",-0.999,[{decimals, 3}]),
     test_fts("1.0",0.999, [{decimals, 2}, compact]),
@@ -184,6 +185,9 @@ t_float_to_string(Config) when is_list(Config) ->
     test_fts("123000000000000000000.0",1.23e20, [{decimals,   10}, compact]),
     test_fts("1.2300000000e+20",1.23e20, [{scientific, 10}, compact]),
     test_fts("1.23000000000000000000e+20",1.23e20, []),
+
+    fts_rand_float_decimals(1000),
+
     ok.
 
 test_fts(Expect, Float) ->
@@ -197,6 +201,49 @@ test_fts(Expect, Float, Args) ->
     BinExpect = float_to_binary(Float,Args).
 
 
+rand_float_reasonable() ->
+    F = rand_float(),
+    case abs(F) > 1.0e238 of
+        true -> rand_float_reasonable();
+        false -> F
+    end.
+
+fts_rand_float_decimals(0) -> ok;
+fts_rand_float_decimals(N) ->
+    [begin
+         F0 = rand_float_reasonable(),
+         L0 = float_to_list(F0, [{decimals, D}]),
+         L1 = case D of
+                  0 -> L0 ++ ".0";
+                  _ -> L0
+              end,
+         F1 = list_to_float(L1),
+         Diff = abs(F0-F1),
+         MaxDiff = max_diff_decimals(F0, D),
+         ok = case Diff =< MaxDiff of
+                  true -> ok;
+                  false ->
+                      io:format("F0 = ~w ~w\n",  [F0, <<F0/float>>]),
+                      io:format("L1 = ~s\n",  [L1]),
+                      io:format("F1 = ~w ~w\n",  [F1, <<F1/float>>]),
+                      io:format("Diff = ~w, MaxDiff = ~w\n", [Diff, MaxDiff]),
+                      error
+              end
+     end
+     || D <- lists:seq(0,15)],
+
+    fts_rand_float_decimals(N-1).
+
+max_diff_decimals(F, D) ->
+    IntBits = floor(math:log2(abs(F))) + 1,
+    FracBits = (52 - IntBits),
+    Log10_2 = 0.3010299956639812,  % math:log10(2)
+    MaxDec = floor(FracBits * Log10_2),
+
+    Resolution = math:pow(2, IntBits - 53),
+
+    (math:pow(10, -min(D,MaxDec)) / 2) + Resolution.
+
 %% Tests list_to_float/1.
 
 t_string_to_float_safe(Config) when is_list(Config) ->
@@ -331,18 +378,26 @@ t_trunc_and_friends(_Config) ->
     -18446744073709551616 = trunc_and_friends(-float(1 bsl 64)),
 
     %% Random.
+    rand_seed(),
     t_trunc_and_friends_rand(100),
     ok.
 
+rand_seed() ->
+    rand:seed(exrop),
+    io:format("\n*** rand:export_seed() = ~w\n\n", [rand:export_seed()]),
+    ok.
+
+rand_float() ->
+    F0 = rand:uniform() * math:pow(10, 50*rand:normal()),
+    case rand:uniform() of
+        U when U < 0.5 -> -F0;
+        _ -> F0
+    end.
+
 t_trunc_and_friends_rand(0) ->
     ok;
 t_trunc_and_friends_rand(N) ->
-    F0 = rand:uniform() * math:pow(10, 50*rand:normal()),
-    F = case rand:uniform() of
-	    U when U < 0.5 -> -F0;
-	    _ -> F0
-	end,
-    _ = trunc_and_friends(F),
+    _ = trunc_and_friends(rand_float()),
     t_trunc_and_friends_rand(N-1).
 
 trunc_and_friends(F) ->