erts: Round bignum to float conversion correctly

erts/emulator/beam/big.c

@@ -1936,20 +1936,118 @@ erts_uint64_array_to_big(Uint **hpp, int neg, int len, Uint64 *array)
 int
 big_to_double(Eterm x, double* resp)
 {
-    double d = 0.0;
     Eterm* xp = big_val(x);
     dsize_t xl = BIG_SIZE(xp);
-    ErtsDigit* s = BIG_V(xp) + xl;
+    ErtsDigit* v = BIG_V(xp);
     short xsgn = BIG_SIGN(xp);
-    double dbase = ((double)(D_MASK)+1);
+    ErtsDigit msd;
+    Uint64 mant;
+    int bitlen, guard, msd_bits, half_bit, exp;
+    ErtsDigit lesser_bits;
+    dsize_t i;
+    double d;
 
-    while (xl--) {
-        d = d * dbase + *--s;
+    ASSERT(xl > 0);
+    msd = v[xl-1];
+    ASSERT(msd != 0);
 
-        if (!erts_isfinite(d)) {
-            return -1;
+    /* Bit length of the most significant digit, and of the whole value. */
+    msd_bits = erts_fit_in_bits_uint(msd);
+    bitlen = (xl-1) * D_EXP + msd_bits;
+
+#if D_EXP == 64
+    ERTS_CT_ASSERT(SMALL_BITS > 53);
+    ASSERT(bitlen > 53);
+#elif D_EXP == 32
+    if (bitlen <= 53) {
+        /*
+         * The value fits the double mantissa exactly, so accumulating
+         * digit by digit cannot round.
+         */
+        ASSERT(xl == 1 || xl == 2);
+
+        d = (double) v[0];
+        if (xl == 2) {
+            const double dbase = ((double)(D_MASK)+1);
+            d += ((double) v[1]) * dbase;
+        }
+        *resp = xsgn ? -d : d;
+        return 0;
+    }
+#endif
+
+    /*
+     * More than 1024 bits is at least 2^1024, above the largest finite double
+     * (2^1024 - 2^971). Reject it here rather than in the loop below, which
+     * may visit every digit to determine rounding. Exactly 1024 bits can still
+     * be finite, so it takes the rounding path, where the erts_isfinite()
+     * check catches a mantissa that carries up to 2^1024.
+     */
+    if (bitlen > 1024) {
+        return -1;
+    }
+
+    /*
+     * Wider than the mantissa, so the result must be rounded. Accumulating
+     * `d = d * base + digit` per digit rounds once per digit and compounds
+     * the error, which can land on the wrong side of the true value; IEEE 754
+     * requires the nearest representable double, ties to even.
+     *
+     * First take the top 54 bits (53 of mantissa plus one half bit).
+     *
+     * Then visit as few lower bignum words as possible to determine rounding.
+     * Round up if the half bit is set and either the mantissa is odd
+     * or some lesser bits are set.
+     */
+    guard = bitlen - 54;
+    mant = 0;
+
+    for (i = xl-1; ; i--) {
+        ErtsDigit dig = v[i];
+        Uint lsb = i * D_EXP;   /* bit position of this digit's LSB */
+
+        if (lsb > guard) {
+            /* Entirely above the cut; every set bit lands within 54 bits. */
+            mant |= (Uint64)dig << (lsb - guard);
+        } else {
+            /* Lowest part of mantissa plus maybe some lesser bits */
+            Uint k = guard - lsb;
+
+            lesser_bits = (dig & (((ErtsDigit)1 << k) - 1));
+            mant |= (Uint64)(dig >> k);
+            break;
+        }
+    }
+
+    half_bit = (int)(mant & 1);
+    mant >>= 1;                      /* 53 significant bits remain */
+    exp = (int)(guard + 1);
+
+    /* Round to nearest, ties to even. */
+    if (half_bit) {
+        if (!(mant & 1)) {
+            while (!lesser_bits) {
+                if (i == 0) {
+                    /* Exactly even and a half, round down */
+                    goto mant_exp_done;
+                }
+                lesser_bits = v[--i];
+            }
+        }
+
+        /* Round up */
+        mant++;
+        if (mant == ((Uint64)1 << 53)) {
+            mant >>= 1;              /* carried out of the mantissa */
+            exp++;
         }
     }
+mant_exp_done:
+
+    d = ldexp((double)mant, exp);
+    if (!erts_isfinite(d)) {
+        return -1;
+    }
 
     *resp = xsgn ? -d : d;
     return 0;

erts/emulator/test/big_SUITE.erl

@@ -24,7 +24,7 @@
 
 -export([t_div/1, eq_28/1, eq_32/1, eq_big/1, eq_math/1, eq_big_mul_div/1,
          big_literals/1, borders/1, negative/1, karatsuba/1,
-         big_float_1/1, big_float_2/1,
+         big_float_1/1, big_float_2/1, big_float_3/1,
          bxor_2pow/1, band_2pow/1,
          shift_limit_1/1, powmod/1, system_limit/1, toobig/1, otp_6692/1,
          properties/1]).
@@ -51,7 +51,7 @@ all() ->
      properties].
 
 groups() -> 
-    [{big_float, [], [big_float_1, big_float_2]}].
+    [{big_float, [], [big_float_1, big_float_2, big_float_3]}].
 
 %%
 %% Syntax of data files:
@@ -338,6 +338,68 @@ big_float_2(Config) when is_list(Config) ->
     {'EXIT', _} = (catch 4/(2*I)),
     ok.
 
+%% Converting a bignum to a float must give the nearest representable
+%% double, ties to even. Accumulating digit by digit rounds once per digit
+%% and compounds the error, which lands on the wrong side of the true value
+%% for some values wider than one digit.
+big_float_3(Config) when is_list(Config) ->
+    rand_seed(),
+    %% Each of these converted to the second-nearest double when the
+    %% conversion rounded per digit.
+    [begin
+         Nearest = correctly_rounded(I),
+         Nearest = float(I),
+         Nearest = 1.0 * I,
+         NegNearest = -Nearest,
+         NegNearest = float(-I)
+     end
+     || I <- [428654966685883400000,
+              38409289721754710000,
+              34784104853086640000,
+              385269108828434300000,
+              96874578115970900000,
+              252558769001389900000,
+              26465126867694860000]],
+
+    %% Widths on both sides of the single-digit boundary, where the
+    %% per-digit accumulation starts to compound.
+    for(50, 300,
+        fun(Bits) ->
+                for(1, 200,
+                    fun(_) ->
+                            I = rand:uniform(1 bsl Bits),
+                            Nearest = correctly_rounded(I),
+                            Nearest = float(I)
+                    end)
+        end),
+
+    %% 2-pows and neighbours
+    [begin
+         I = (1 bsl E) + Diff,
+         Nearest = correctly_rounded(I),
+         Nearest = float(I)
+     end
+     || E <- lists:seq(0, 1023), Diff <- lists:seq(-2,2)],
+
+    %% Mantissa rounding edge cases
+    [begin
+         Mant = (1 bsl 52) + Odd,
+         I = (Mant bsl Exp) + (Half bsl (Exp-1)) + (1 bsl Low),
+         Nearest = correctly_rounded(I),
+         Nearest = float(I)
+     end
+     || Exp <- lists:seq(1, 1023-53),
+        Odd <- [0,1],
+        Half <- [1],
+        Low <- [-1 | lists:seq(0, Exp-2, 8)]],
+
+    ok.
+
+%% The platform's decimal parser is correctly rounded, so it serves as the
+%% oracle for what float/1 must return.
+correctly_rounded(I) ->
+    binary_to_float(iolist_to_binary([integer_to_list(I), ".0"])).
+
 %% OTP-3256
 shift_limit_1(Config) when is_list(Config) ->
     case catch (id(1) bsl 100000000) of