Starting work on a possible integration of dragonbox.

include/fast_float/dragonbox.h

@@ -0,0 +1,208 @@
+#ifndef DRAGONBOX_H
+#define DRAGONBOX_H
+
+#include <cstdint>
+
+/**
+ * Our objective is to integrate a fast binary to decimal converter which relies
+ * on our precomputed constants.
+ */
+
+namespace fast_float {
+// credit: https://github.com/jk-jeon/dragonbox
+// credit: https://github.com/fmtlib/fmt/blob/012cc709d0abb0a963591413a746b988168e5e3a/include/fmt/format-inl.h
+namespace dragonbox {
+
+
+
+// The main algorithm for shorter interval case
+template <class T>
+FMT_INLINE decimal_fp<T> shorter_interval_case(int exponent) FMT_NOEXCEPT {
+  decimal_fp<T> ret_value;
+  // Compute k and beta
+  const int minus_k = floor_log10_pow2_minus_log10_4_over_3(exponent);
+  const int beta_minus_1 = exponent + floor_log2_pow10(-minus_k);
+
+  // Compute xi and zi
+  using cache_entry_type = typename cache_accessor<T>::cache_entry_type;
+  const cache_entry_type cache = cache_accessor<T>::get_cached_power(-minus_k);
+
+  auto xi = cache_accessor<T>::compute_left_endpoint_for_shorter_interval_case(
+      cache, beta_minus_1);
+  auto zi = cache_accessor<T>::compute_right_endpoint_for_shorter_interval_case(
+      cache, beta_minus_1);
+
+  // If the left endpoint is not an integer, increase it
+  if (!is_left_endpoint_integer_shorter_interval<T>(exponent)) ++xi;
+
+  // Try bigger divisor
+  ret_value.significand = zi / 10;
+
+  // If succeed, remove trailing zeros if necessary and return
+  if (ret_value.significand * 10 >= xi) {
+    ret_value.exponent = minus_k + 1;
+    ret_value.exponent += remove_trailing_zeros(ret_value.significand);
+    return ret_value;
+  }
+
+  // Otherwise, compute the round-up of y
+  ret_value.significand =
+      cache_accessor<T>::compute_round_up_for_shorter_interval_case(
+          cache, beta_minus_1);
+  ret_value.exponent = minus_k;
+
+  // When tie occurs, choose one of them according to the rule
+  if (exponent >= float_info<T>::shorter_interval_tie_lower_threshold &&
+      exponent <= float_info<T>::shorter_interval_tie_upper_threshold) {
+    ret_value.significand = ret_value.significand % 2 == 0
+                                ? ret_value.significand
+                                : ret_value.significand - 1;
+  } else if (ret_value.significand < xi) {
+    ++ret_value.significand;
+  }
+  return ret_value;
+}
+
+
+
+/**
+ * This is the main function.
+ */
+template <typename T> decimal_fp<T> to_decimal(T x) noexcept {
+  // Step 1: integer promotion & Schubfach multiplier calculation.
+
+  using carrier_uint = typename float_info<T>::carrier_uint;
+  using cache_entry_type = typename cache_accessor<T>::cache_entry_type;
+  auto br = bit_cast<carrier_uint>(x);
+
+  // Extract significand bits and exponent bits.
+  const carrier_uint significand_mask =
+      (static_cast<carrier_uint>(1) << float_info<T>::significand_bits) - 1;
+  carrier_uint significand = (br & significand_mask);
+  int exponent = static_cast<int>((br & exponent_mask<T>()) >>
+                                  float_info<T>::significand_bits);
+
+  if (exponent != 0) {  // Check if normal.
+    exponent += float_info<T>::exponent_bias - float_info<T>::significand_bits;
+
+    // Shorter interval case; proceed like Schubfach.
+    if (significand == 0) return shorter_interval_case<T>(exponent);
+
+    significand |=
+        (static_cast<carrier_uint>(1) << float_info<T>::significand_bits);
+  } else {
+    // Subnormal case; the interval is always regular.
+    if (significand == 0) return {0, 0};
+    exponent = float_info<T>::min_exponent - float_info<T>::significand_bits;
+  }
+
+  const bool include_left_endpoint = (significand % 2 == 0);
+  const bool include_right_endpoint = include_left_endpoint;
+
+  // Compute k and beta.
+  const int minus_k = floor_log10_pow2(exponent) - float_info<T>::kappa;
+  const cache_entry_type cache = cache_accessor<T>::get_cached_power(-minus_k);
+  const int beta_minus_1 = exponent + floor_log2_pow10(-minus_k);
+
+  // Compute zi and deltai
+  // 10^kappa <= deltai < 10^(kappa + 1)
+  const uint32_t deltai = cache_accessor<T>::compute_delta(cache, beta_minus_1);
+  const carrier_uint two_fc = significand << 1;
+  const carrier_uint two_fr = two_fc | 1;
+  const carrier_uint zi =
+      cache_accessor<T>::compute_mul(two_fr << beta_minus_1, cache);
+
+  // Step 2: Try larger divisor; remove trailing zeros if necessary
+
+  // Using an upper bound on zi, we might be able to optimize the division
+  // better than the compiler; we are computing zi / big_divisor here
+  decimal_fp<T> ret_value;
+  ret_value.significand = divide_by_10_to_kappa_plus_1(zi);
+  uint32_t r = static_cast<uint32_t>(zi - float_info<T>::big_divisor *
+                                              ret_value.significand);
+
+  if (r > deltai) {
+    goto small_divisor_case_label;
+  } else if (r < deltai) {
+    // Exclude the right endpoint if necessary
+    if (r == 0 && !include_right_endpoint &&
+        is_endpoint_integer<T>(two_fr, exponent, minus_k)) {
+      --ret_value.significand;
+      r = float_info<T>::big_divisor;
+      goto small_divisor_case_label;
+    }
+  } else {
+    // r == deltai; compare fractional parts
+    // Check conditions in the order different from the paper
+    // to take advantage of short-circuiting
+    const carrier_uint two_fl = two_fc - 1;
+    if ((!include_left_endpoint ||
+         !is_endpoint_integer<T>(two_fl, exponent, minus_k)) &&
+        !cache_accessor<T>::compute_mul_parity(two_fl, cache, beta_minus_1)) {
+      goto small_divisor_case_label;
+    }
+  }
+  ret_value.exponent = minus_k + float_info<T>::kappa + 1;
+
+  // We may need to remove trailing zeros
+  ret_value.exponent += remove_trailing_zeros(ret_value.significand);
+  return ret_value;
+
+  // Step 3: Find the significand with the smaller divisor
+
+small_divisor_case_label:
+  ret_value.significand *= 10;
+  ret_value.exponent = minus_k + float_info<T>::kappa;
+
+  const uint32_t mask = (1u << float_info<T>::kappa) - 1;
+  auto dist = r - (deltai / 2) + (float_info<T>::small_divisor / 2);
+
+  // Is dist divisible by 2^kappa?
+  if ((dist & mask) == 0) {
+    const bool approx_y_parity =
+        ((dist ^ (float_info<T>::small_divisor / 2)) & 1) != 0;
+    dist >>= float_info<T>::kappa;
+
+    // Is dist divisible by 5^kappa?
+    if (check_divisibility_and_divide_by_pow5<float_info<T>::kappa>(dist)) {
+      ret_value.significand += dist;
+
+      // Check z^(f) >= epsilon^(f)
+      // We have either yi == zi - epsiloni or yi == (zi - epsiloni) - 1,
+      // where yi == zi - epsiloni if and only if z^(f) >= epsilon^(f)
+      // Since there are only 2 possibilities, we only need to care about the
+      // parity. Also, zi and r should have the same parity since the divisor
+      // is an even number
+      if (cache_accessor<T>::compute_mul_parity(two_fc, cache, beta_minus_1) !=
+          approx_y_parity) {
+        --ret_value.significand;
+      } else {
+        // If z^(f) >= epsilon^(f), we might have a tie
+        // when z^(f) == epsilon^(f), or equivalently, when y is an integer
+        if (is_center_integer<T>(two_fc, exponent, minus_k)) {
+          ret_value.significand = ret_value.significand % 2 == 0
+                                      ? ret_value.significand
+                                      : ret_value.significand - 1;
+        }
+      }
+    }
+    // Is dist not divisible by 5^kappa?
+    else {
+      ret_value.significand += dist;
+    }
+  }
+  // Is dist not divisible by 2^kappa?
+  else {
+    // Since we know dist is small, we might be able to optimize the division
+    // better than the compiler; we are computing dist / small_divisor here
+    ret_value.significand +=
+        small_division_by_pow10<float_info<T>::kappa>(dist);
+  }
+  return ret_value;
+}
+
+}
+} // fast_float
+
+
+#endif