From 4d6f28499d4514eed337be1b43233a89455a9eb9 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Cl=C3=A9ment=20Joubert?= Date: Mon, 26 Oct 2026 19:56:52 +0300 Subject: [PATCH 1/3] Faster Float to BigDecimal conversion for doubles with a short decimal form Doubles reading back from m / 10**k with m < 2**51 and k <= 22 (prices, rates, most decimal literals) get their shortest digits without dtoa, which allocates and frees several Bigints per call since #528. Co-Authored-By: Claude Opus 5.5 --- ext/bigdecimal/bigdecimal.c | 73 +++++++++++++++++++++---- test/bigdecimal/test_bigdecimal_util.rb | 19 +++++++ 2 files changed, 82 insertions(+), 10 deletions(-) diff --git a/ext/bigdecimal/bigdecimal.c b/ext/bigdecimal/bigdecimal.c index e934e340..a032d0f5 100644 --- a/ext/bigdecimal/bigdecimal.c +++ b/ext/bigdecimal/bigdecimal.c @@ -2662,6 +2662,53 @@ rb_inum_convert_to_BigDecimal(VALUE val) } } +/* + * Shortest round-trip digits of a finite nonzero double, as BigDecimal_dtoa + * mode 0 returns them, for doubles that read back from m / 10**k with + * m < 2**51 and 0 <= k <= 22: prices, rates and most decimal literals. + * + * At such a scale, the doubles reading back to d span less than 0.5 units, + * so at most one integer m reads back to d. The smallest k with such an m + * gives the shortest digits, and IEEE division of the two exact operands + * rounds m / 10**k exactly like strtod would. + * + * Returns the number of digits written to buf, or 0 when dtoa is needed. + */ +static int +float_shortest_digits_fast(double d, char *buf, int *decpt, int *negative_p) +{ +#if defined(FLT_EVAL_METHOD) && FLT_EVAL_METHOD == 0 + static const double powers_of_ten[] = { + 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9, 1e10, 1e11, + 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19, 1e20, 1e21, 1e22 + }; + const double two_pow_51 = 2251799813685248.0; + const double abs_d = fabs(d); + + for (int k = 0; k < (int)numberof(powers_of_ten); k++) { + const double scaled = abs_d * powers_of_ten[k]; + if (!(scaled < two_pow_51)) return 0; + + const double m = round(scaled); + if (m / powers_of_ten[k] != abs_d) continue; + + char digits[BIGDECIMAL_DOUBLE_FIGURES + 1]; + uint64_t n = (uint64_t)m; + int len = 0; + for (; n > 0; n /= 10) digits[len++] = '0' + (char)(n % 10); + + int ntz = 0; + while (digits[ntz] == '0') ntz++; + for (int i = 0; i < len - ntz; i++) buf[i] = digits[len - 1 - i]; + + *decpt = len - k; + *negative_p = signbit(d) ? 1 : 0; + return len - ntz; + } +#endif + return 0; +} + static VALUE rb_float_convert_to_BigDecimal(VALUE val, size_t digs, int raise_exception) { @@ -2704,16 +2751,22 @@ rb_float_convert_to_BigDecimal(VALUE val, size_t digs, int raise_exception) /* Use the same logic in flo_to_s to convert a float to a decimal string */ char buf[BIGDECIMAL_DOUBLE_FIGURES + BASE_FIG + 2 + 1]; /* sizeof(buf) == 28 in the typical case */ int decpt, negative_p; - char *e; - const int mode = digs == 0 ? 0 : 2; - char *p = BigDecimal_dtoa(d, mode, (int)digs, &decpt, &negative_p, &e); - int len10 = (int)(e - p); - if (len10 > BIGDECIMAL_DOUBLE_FIGURES) { - /* TODO: Presumably, rounding should be done here. */ - len10 = BIGDECIMAL_DOUBLE_FIGURES; - } - memcpy(buf, p, len10); - free(p); + int len10 = 0; + if (digs == 0) { + len10 = float_shortest_digits_fast(d, buf, &decpt, &negative_p); + } + if (len10 == 0) { + char *e; + const int mode = digs == 0 ? 0 : 2; + char *p = BigDecimal_dtoa(d, mode, (int)digs, &decpt, &negative_p, &e); + len10 = (int)(e - p); + if (len10 > BIGDECIMAL_DOUBLE_FIGURES) { + /* TODO: Presumably, rounding should be done here. */ + len10 = BIGDECIMAL_DOUBLE_FIGURES; + } + memcpy(buf, p, len10); + free(p); + } VALUE inum; size_t RB_UNUSED_VAR(prec) = 0; diff --git a/test/bigdecimal/test_bigdecimal_util.rb b/test/bigdecimal/test_bigdecimal_util.rb index 67693935..5e1d951b 100644 --- a/test/bigdecimal/test_bigdecimal_util.rb +++ b/test/bigdecimal/test_bigdecimal_util.rb @@ -40,6 +40,25 @@ def test_Float_to_d_without_precision assert_equal(BigDecimal("999_999.9999"), 999_999.9999.to_d) end + def test_Float_to_d_shortest_representation + floats = [ + 0.1, 0.2, 129.9, 108.25, 19.99, 1.2, 0.055, 5e-324, 1e-22, 1.5e-22, 1e-23, 1e22, 1e23, + 2.0**-30, 2.0**50, 2.0**51 - 1, 2.0**51 + 2, 2.0**52 + 1, 2.0**53 - 1, + (2.0**51 - 1) / 10, 225179981368524.7, 9007199254740.993, + 0.1.prev_float, 0.1.next_float, 129.9.prev_float, 1.0.prev_float, 1.0.next_float + ] + random = Random.new(42) + floats += 1000.times.map { random.rand(1_000_000) / 10.0**random.rand(0..22) } + floats += 1000.times.map { random.bytes(8).unpack1("D") }.select(&:finite?) + floats.each do |float| + [float, -float].each do |f| + significand = f.to_s.split("e").first.delete("^0-9").sub(/\A0+/, "").sub(/0+\z/, "") + next if significand.size > BigDecimal.double_fig # Float#to_d keeps double_fig digits only + assert_equal(BigDecimal(f.to_s), f.to_d, "#{f.inspect}.to_d") + end + end + end + def test_Float_to_d_with_precision digits = 5 delta = 1.0/10**(digits) From 123cd59af5323a27e902db0feba27cfe290ed29a Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Cl=C3=A9ment=20Joubert?= Date: Tue, 27 Oct 2026 09:42:11 +0300 Subject: [PATCH 2/3] Skip the Float#to_d fast path under -ffast-math and below 1e-22 -ffast-math keeps FLT_EVAL_METHOD at 0 but lets the compiler turn m / 10**k into m * 10**-k, which breaks the exact division the fast path relies on (1.0.prev_float.to_d would give 1 on x86_64 clang -O3). No double below 1e-22 can read back from m / 10**k with k <= 22, so return early instead of running 23 multiplications that may underflow. Co-Authored-By: Claude Opus 5.5 --- ext/bigdecimal/bigdecimal.c | 5 +++-- 1 file changed, 3 insertions(+), 2 deletions(-) diff --git a/ext/bigdecimal/bigdecimal.c b/ext/bigdecimal/bigdecimal.c index a032d0f5..056000ee 100644 --- a/ext/bigdecimal/bigdecimal.c +++ b/ext/bigdecimal/bigdecimal.c @@ -2677,15 +2677,16 @@ rb_inum_convert_to_BigDecimal(VALUE val) static int float_shortest_digits_fast(double d, char *buf, int *decpt, int *negative_p) { -#if defined(FLT_EVAL_METHOD) && FLT_EVAL_METHOD == 0 +#if defined(FLT_EVAL_METHOD) && FLT_EVAL_METHOD == 0 && !defined(__FAST_MATH__) static const double powers_of_ten[] = { 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9, 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19, 1e20, 1e21, 1e22 }; const double two_pow_51 = 2251799813685248.0; const double abs_d = fabs(d); + if (abs_d < 1e-22) return 0; - for (int k = 0; k < (int)numberof(powers_of_ten); k++) { + for (int k = 0; k < numberof(powers_of_ten); k++) { const double scaled = abs_d * powers_of_ten[k]; if (!(scaled < two_pow_51)) return 0; From 659c937750e8b1f2e00125cb9222ada437750e42 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Cl=C3=A9ment=20Joubert?= Date: Tue, 27 Oct 2026 12:33:36 +0300 Subject: [PATCH 3/3] Verify Float#to_d fast path candidates in exact integer arithmetic The fast path accepted m / 10**k when the FP division gave back d, which only proves round-trip in round-to-nearest without reciprocal math. Under FE_DOWNWARD 0.1.prev_float.to_d gave 0.1, and an x86_64 -O3 -ffast-math build turned 1.0.prev_float.to_d into 1. The division stays as a cheap filter, but a candidate is now accepted only when |m * 2**q - 2 * mantissa * 5**k| < 5**k holds in 128-bit integers, so the FPU rounding mode and fast-math can at worst send a double to dtoa. This replaces the FLT_EVAL_METHOD and __FAST_MATH__ gates with HAVE_UINT128_T. Co-Authored-By: Claude Opus 5.5 --- ext/bigdecimal/bigdecimal.c | 54 +++++++++++++++++++++++++++++++------ 1 file changed, 46 insertions(+), 8 deletions(-) diff --git a/ext/bigdecimal/bigdecimal.c b/ext/bigdecimal/bigdecimal.c index 056000ee..1e4bfe6f 100644 --- a/ext/bigdecimal/bigdecimal.c +++ b/ext/bigdecimal/bigdecimal.c @@ -2662,22 +2662,54 @@ rb_inum_convert_to_BigDecimal(VALUE val) } } +#ifdef HAVE_UINT128_T +/* + * Whether m / 10**k rounds to nearest to mantissa * 2**exponent, a normal + * double. The comparison is exact, so it holds whatever the FPU rounding + * mode or -ffast-math: with q = 1 - exponent - k, the condition + * |m / 10**k - d| < ulp(d) / 2 is |m * 2**q - 2 * mantissa * 5**k| < 5**k, + * and the gap below a power of two is half as wide. Ties fail; they cannot + * occur for m < 2**52, k <= 22. + */ +static int +float_reads_back_from_decimal(uint64_t mantissa, int exponent, uint64_t m, int k) +{ + static const uint64_t powers_of_five[] = { + 1ULL, 5ULL, 25ULL, 125ULL, + 625ULL, 3125ULL, 15625ULL, 78125ULL, + 390625ULL, 1953125ULL, 9765625ULL, 48828125ULL, + 244140625ULL, 1220703125ULL, 6103515625ULL, 30517578125ULL, + 152587890625ULL, 762939453125ULL, 3814697265625ULL, 19073486328125ULL, + 95367431640625ULL, 476837158203125ULL, 2384185791015625ULL + }; + const int q = 1 - exponent - k; + if (q < 0 || q > 127 || ((uint128_t)m >> (127 - q)) != 0) return 0; + + const uint128_t decimal = (uint128_t)m << q; + const uint128_t binary = (uint128_t)(2 * mantissa) * powers_of_five[k]; + const uint128_t half_ulp = powers_of_five[k]; + if (decimal >= binary) return decimal - binary < half_ulp; + if (mantissa == 1ULL << (DBL_MANT_DIG - 1)) return 2 * (binary - decimal) < half_ulp; + return binary - decimal < half_ulp; +} +#endif + /* * Shortest round-trip digits of a finite nonzero double, as BigDecimal_dtoa * mode 0 returns them, for doubles that read back from m / 10**k with - * m < 2**51 and 0 <= k <= 22: prices, rates and most decimal literals. + * m <= 2**51 and 0 <= k <= 22: prices, rates and most decimal literals. * * At such a scale, the doubles reading back to d span less than 0.5 units, - * so at most one integer m reads back to d. The smallest k with such an m - * gives the shortest digits, and IEEE division of the two exact operands - * rounds m / 10**k exactly like strtod would. + * so at most one integer m reads back to d, and the nearest integer to + * d * 10**k is it whenever it exists. Every k that finds it gives the same + * digits once trailing zeros are stripped: the shortest ones. * * Returns the number of digits written to buf, or 0 when dtoa is needed. */ static int float_shortest_digits_fast(double d, char *buf, int *decpt, int *negative_p) { -#if defined(FLT_EVAL_METHOD) && FLT_EVAL_METHOD == 0 && !defined(__FAST_MATH__) +#ifdef HAVE_UINT128_T static const double powers_of_ten[] = { 1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9, 1e10, 1e11, 1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19, 1e20, 1e21, 1e22 @@ -2686,15 +2718,21 @@ float_shortest_digits_fast(double d, char *buf, int *decpt, int *negative_p) const double abs_d = fabs(d); if (abs_d < 1e-22) return 0; + int exponent; + const uint64_t mantissa = (uint64_t)ldexp(frexp(abs_d, &exponent), DBL_MANT_DIG); + exponent -= DBL_MANT_DIG; + for (int k = 0; k < numberof(powers_of_ten); k++) { const double scaled = abs_d * powers_of_ten[k]; if (!(scaled < two_pow_51)) return 0; - const double m = round(scaled); - if (m / powers_of_ten[k] != abs_d) continue; + const double rounded = round(scaled); + if (rounded / powers_of_ten[k] != abs_d) continue; /* cheap filter, the exact check decides */ + const uint64_t m = (uint64_t)rounded; + if (!float_reads_back_from_decimal(mantissa, exponent, m, k)) continue; char digits[BIGDECIMAL_DOUBLE_FIGURES + 1]; - uint64_t n = (uint64_t)m; + uint64_t n = m; int len = 0; for (; n > 0; n /= 10) digits[len++] = '0' + (char)(n % 10);