From aae60d04ee2988fc349b3e081e76289b21f7e876 Mon Sep 17 00:00:00 2001 From: Stanislav Yaglo Date: Fri, 24 Jul 2026 23:15:27 +0100 Subject: [PATCH] erts: Round bignum to float conversion correctly MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit big_to_double/2 accumulated the result one digit at a time: while (xl--) { d = d * dbase + *--s; } Each iteration rounds, so the errors compound and the result can be the second-nearest double rather than the nearest. IEEE 754 requires the nearest representable value, ties to even. A value that fits in a single digit rounds only once and was already correct, so the defect appears from two digits up: 16% of 65-bit integers converted to the wrong double, and never to the nearer one. For example, float(428654966685883400000) returned 4.2865496668588343e20 where the nearest double is 4.286549666858834e20. The affected value is also what binary_to_float/1 returns for the same digits, since that routes through the platform's correctly rounded decimal parser. Compute the bit length instead, and for values wider than the 53-bit mantissa take the top 54 bits — 53 of mantissa plus one round bit. Rounding is then determined by reading as few of the lower digits as possible. Values of 53 bits or fewer keep the digit-by-digit accumulation, which cannot round there. Integers wider than 1024 bits return -1 up front: they are at least 2^1024 and cannot be finite. The previous loop stopped as soon as the accumulator went infinite, and without this check a conversion would become O(size) for bignums that reach tens of thousands of digits. Exactly 1024 bits can still be finite and takes the rounding path, where the overflow check after ldexp() catches a mantissa that carries up to 2^1024. The defect surfaced while implementing RFC 8785 (JSON Canonicalization Scheme), which serializes numbers as ECMAScript does and therefore has to decide whether an integer's decimal digits survive a round trip through a double. Comparing float/1 against binary_to_float/1 over random integers disagreed often enough to be reproducible, and the disagreement was always in binary_to_float/1's favour. big_float_3 covers the reproducers, random values from 50 to 300 bits, and every power of two up to 2^1023; it fails on the previous implementation. Co-authored-by: Sverker Eriksson --- erts/emulator/beam/big.c | 112 +++++++++++++++++++++++++++++-- erts/emulator/test/big_SUITE.erl | 66 +++++++++++++++++- 2 files changed, 169 insertions(+), 9 deletions(-) diff --git a/erts/emulator/beam/big.c b/erts/emulator/beam/big.c index e3a12e324c5e..0a98d7716e58 100644 --- a/erts/emulator/beam/big.c +++ b/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; diff --git a/erts/emulator/test/big_SUITE.erl b/erts/emulator/test/big_SUITE.erl index 5e121e30c570..aecdce0420b9 100644 --- a/erts/emulator/test/big_SUITE.erl +++ b/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