erts: Round bignum to float conversion correctly

erts/emulator/beam/big.c

@@ -2411,19 +2411,99 @@ 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 bitlen, guard, mant;
+    int msd_bits, sticky, round_bit, exp;
+    dsize_t i;
+    double d;
+
+    ASSERT(xl > 0);
+    msd = v[xl-1];
+    ASSERT(msd != 0);
+
+    /* Bit length of the most significant digit, and of the whole value. */
+    msd_bits = erts_fit_in_bits_uint(msd);
+    bitlen = (Uint64)(xl-1) * D_EXP + msd_bits;
 
-    while (xl--) {
-        d = d * dbase + *--s;
+    if (bitlen <= 53) {
+        /*
+         * The value fits the double mantissa exactly, so accumulating
+         * digit by digit cannot round.
+         */
+        double dbase = ((double)(D_MASK)+1);
 
-        if (!erts_isfinite(d)) {
-            return -1;
+        d = 0.0;
+        for (i = xl; i-- > 0; ) {
+            d = d * dbase + v[i];
         }
+        *resp = xsgn ? -d : d;
+        return 0;
+    }
+
+    /*
+     * 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
+     * visits every digit to collect the sticky bit; bignums reach tens of
+     * thousands of digits. 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.
+     *
+     * Take the top 54 bits (53 of mantissa plus one round bit) and record
+     * whether any bit below them is set, then round exactly once.
+     */
+    guard = bitlen - 54;
+    mant = 0;
+    sticky = 0;
+
+    for (i = xl; i-- > 0 && !sticky; ) {
+        ErtsDigit dig = v[i];
+        Uint64 lsb = (Uint64)i * D_EXP;   /* bit position of this digit's LSB */
+
+        if (lsb + D_EXP <= guard) {
+            /* Entirely below the cut. */
+            sticky |= (dig != 0);
+        } else if (lsb >= guard) {
+            /* Entirely above the cut; every set bit lands within 54 bits. */
+            mant |= (Uint64)dig << (lsb - guard);
+        } else {
+            /* Straddles the cut. */
+            Uint64 k = guard - lsb;
+
+            sticky |= (dig & (((ErtsDigit)1 << k) - 1)) != 0;
+            mant |= (Uint64)(dig >> k);
+        }
+    }
+
+    round_bit = (int)(mant & 1);
+    mant >>= 1;                      /* 53 significant bits remain */
+    exp = (int)(guard + 1);
+
+    /* Round to nearest, ties to even. */
+    if (round_bit && (sticky || (mant & 1))) {
+        mant++;
+        if (mant == ((Uint64)1 << 53)) {
+            mant >>= 1;              /* carried out of the mantissa */
+            exp++;
+        }
+    }
+
+    d = ldexp((double)mant, exp);
+    if (!erts_isfinite(d)) {
+        return -1;
     }
 
     *resp = xsgn ? -d : d;

erts/emulator/test/big_SUITE.erl

@@ -27,7 +27,7 @@
 -export([t_div/1, eq_28/1, eq_32/1, eq_big/1, eq_math/1, eq_big_mul_div/1,
          eq_big_rem/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, reductions/1]).
@@ -56,7 +56,7 @@ all() ->
      properties, reductions].
 
 groups() -> 
-    [{big_float, [], [big_float_1, big_float_2]}].
+    [{big_float, [], [big_float_1, big_float_2, big_float_3]}].
 
 %%
 %% Syntax of data files:
@@ -353,6 +353,52 @@ 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) ->
+    %% 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.
+    [begin
+         I = rand:uniform(1 bsl Bits),
+         Nearest = correctly_rounded(I),
+         Nearest = float(I)
+     end
+     || Bits <- lists:seq(50, 300), _ <- lists:seq(1, 2000)],
+
+    %% 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)],
+
+    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