Float#to_d 2x faster for short decimals like prices and rates, skipping dtoa

ext/bigdecimal/bigdecimal.c

@@ -2662,6 +2662,92 @@ 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.
+ *
+ * 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, 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)
+{
+#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
+    };
+    const double two_pow_51 = 2251799813685248.0;
+    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 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 = 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 +2790,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;

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)