Implement handmade FP

include/fmt/format-inl.h

@@ -219,8 +219,9 @@ FMT_FUNC void system_error::init(
   base = std::runtime_error(to_string(buffer));
 }
 
+namespace internal {
 template <typename T>
-int internal::char_traits<char>::format_float(
+int char_traits<char>::format_float(
     char *buffer, std::size_t size, const char *format,
     unsigned width, int precision, T value) {
   if (width == 0) {
@@ -234,7 +235,7 @@ int internal::char_traits<char>::format_float(
 }
 
 template <typename T>
-int internal::char_traits<wchar_t>::format_float(
+int char_traits<wchar_t>::format_float(
     wchar_t *buffer, std::size_t size, const wchar_t *format,
     unsigned width, int precision, T value) {
   if (width == 0) {
@@ -248,7 +249,7 @@ int internal::char_traits<wchar_t>::format_float(
 }
 
 template <typename T>
-const char internal::basic_data<T>::DIGITS[] =
+const char basic_data<T>::DIGITS[] =
     "0001020304050607080910111213141516171819"
     "2021222324252627282930313233343536373839"
     "4041424344454647484950515253545556575859"
@@ -267,18 +268,30 @@ const char internal::basic_data<T>::DIGITS[] =
   factor * 1000000000
 
 template <typename T>
-const uint32_t internal::basic_data<T>::POWERS_OF_10_32[] = {
+const uint32_t basic_data<T>::POWERS_OF_10_32[] = {
   0, FMT_POWERS_OF_10(1)
 };
 
 template <typename T>
-const uint64_t internal::basic_data<T>::POWERS_OF_10_64[] = {
+const uint64_t basic_data<T>::POWERS_OF_10_64[] = {
   0,
   FMT_POWERS_OF_10(1),
   FMT_POWERS_OF_10(1000000000ull),
   10000000000000000000ull
 };
 
+FMT_FUNC fp operator*(fp x, fp y) {
+  // Multiply 32-bit parts of significands.
+  uint64_t mask = (1ULL << 32) - 1;
+  uint64_t a = x.f >> 32, b = x.f & mask;
+  uint64_t c = y.f >> 32, d = y.f & mask;
+  uint64_t ac = a * c, bc = b * c, ad = a * d, bd = b * d;
+  // Compute mid 64-bit of result and round.
+  uint64_t mid = (bd >> 32) + (ad & mask) + (bc & mask) + (1U << 31);
+  return fp(ac + (ad >> 32) + (bc >> 32) + (mid >> 32), x.e + y.e + 64); 
+}
+}  // namespace internal
+
 #if FMT_USE_WINDOWS_H
 
 FMT_FUNC internal::utf8_to_utf16::utf8_to_utf16(string_view s) {

include/fmt/format.h

@@ -245,6 +245,30 @@ inline dummy_int _finite(...) { return dummy_int(); }
 inline dummy_int isnan(...) { return dummy_int(); }
 inline dummy_int _isnan(...) { return dummy_int(); }
 
+// A handmade floating-point number f * pow(2, e).
+struct fp {
+  uint64_t f;
+  int e;
+  fp(uint64_t f, int e): f(f), e(e) {}
+};
+
+// Returns an fp number representing x - y. Result may not be normalized.
+inline fp operator-(fp x, fp y) {
+  FMT_ASSERT(x.f >= y.f && x.e == y.e, "invalid operands");
+  return fp(x.f - y.f, x.e); 
+}
+
+// Computes an fp number r with r.f = x.f * y.f / pow(2, 32) rounded to nearest
+// with half-up tie breaking, r.e = x.e + y.e + 32. Result may not be normalized.
+fp operator*(fp x, fp y);
+
+// Compute k such that its cached power c_k = c_k.f * pow(2, c_k.e) satisfies
+// alpha <= c_k.e + e <= alpha + 3.
+inline int compute_cached_power_index(int e, int alpha) {
+  constexpr double one_over_log2_10 = 0.30102999566398114;  // 1 / log2(10)
+  return std::ceil((alpha - e + 63) * one_over_log2_10);
+}
+
 template <typename Allocator>
 typename Allocator::value_type *allocate(Allocator& alloc, std::size_t n) {
 #if __cplusplus >= 201103L || FMT_MSC_VER >= 1700

test/util-test.cc

@@ -37,6 +37,7 @@ using fmt::basic_format_arg;
 using fmt::internal::basic_buffer;
 using fmt::basic_memory_buffer;
 using fmt::string_view;
+using fmt::internal::fp;
 using fmt::internal::value;
 
 using testing::_;
@@ -869,3 +870,18 @@ TEST(UtilTest, ParseNonnegativeInt) {
         parse_nonnegative_int(s, fmt::internal::error_handler()),
         fmt::format_error, "number is too big");
 }
+
+TEST(UtilTest, FPSubtract) {
+  auto r = fp(123, 1) - fp(102, 1);
+  EXPECT_EQ(r.f, 21u);
+  EXPECT_EQ(r.e, 1);
+}
+
+TEST(UtilTest, FPMultiply) {
+  auto r = fp(123ULL << 32, 4) * fp(56ULL << 32, 7);
+  EXPECT_EQ(r.f, 123u * 56u);
+  EXPECT_EQ(r.e, 4 + 7 + 64);
+  r = fp(123ULL << 32, 4) * fp(567ULL << 31, 8);
+  EXPECT_EQ(r.f, (123 * 567 + 1u) / 2);
+  EXPECT_EQ(r.e, 4 + 8 + 64);
+}