diff --git a/erts/emulator/beam/big.c b/erts/emulator/beam/big.c index 943f0adcc79d..1c8f75ae1ff2 100644 --- a/erts/emulator/beam/big.c +++ b/erts/emulator/beam/big.c @@ -2411,19 +2411,116 @@ 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; + + ASSERT(xl > 0); + msd = v[xl-1]; + ASSERT(msd != 0); - while (xl--) { - d = d * dbase + *--s; + /* 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); - if (!erts_isfinite(d)) { - return -1; + 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; true; 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 */ + goto no_rounding; + } + lesser_bits = v[--i]; + } + } + + mant++; + if (mant == ((Uint64)1 << 53)) { + mant >>= 1; /* carried out of the mantissa */ + exp++; + } + } +no_rounding: + + d = ldexp((double)mant, exp); + if (!erts_isfinite(d)) { + return -1; } *resp = xsgn ? -d : d; diff --git a/erts/emulator/test/big_SUITE.erl b/erts/emulator/test/big_SUITE.erl index 2533f5b69f2f..9bffb470d805 100644 --- a/erts/emulator/test/big_SUITE.erl +++ b/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,56 @@ 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, 2000, + 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)], + + 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