From c41cea8b18232948f8970eef2efe62670c72b396 Mon Sep 17 00:00:00 2001
From: Victor Zverovich <victor.zverovich@gmail.com>
Date: Sat, 5 Oct 2019 11:45:33 -0700
Subject: [PATCH] Initial implementation of square

---
 include/fmt/format-inl.h | 66 ++++++++++++++++++++++++++++++++++------
 test/format-impl-test.cc | 32 +++++++++++++++++++
 2 files changed, 88 insertions(+), 10 deletions(-)

diff --git a/include/fmt/format-inl.h b/include/fmt/format-inl.h
index 62de4bff..ddf66646 100644
--- a/include/fmt/format-inl.h
+++ b/include/fmt/format-inl.h
@@ -405,7 +405,7 @@ class fp {
     }
     // Subtract 1 to account for hidden bit.
     const auto offset =
-      fp::significand_size - fp::double_significand_size - SHIFT - 1;
+        fp::significand_size - fp::double_significand_size - SHIFT - 1;
     value.f <<= offset;
     value.e -= offset;
     return value;
@@ -468,6 +468,27 @@ FMT_FUNC fp get_cached_power(int min_exponent, int& pow10_exponent) {
   return fp(data::pow10_significands[index], data::pow10_exponents[index]);
 }
 
+// A simple accumulator to hold the sums of terms in bigint::square if uint128_t
+// is not available.
+struct accumulator {
+  uint64_t lower;
+  uint64_t upper;
+
+  accumulator() : lower(0), upper(0) {}
+  explicit operator uint32_t() const { return static_cast<uint32_t>(lower); }
+
+  void operator+=(uint64_t n) {
+    lower += n;
+    if (lower < n) ++upper;
+  }
+  void operator>>=(int shift) {
+    assert(shift == 32);
+    (void)shift;
+    lower = (upper << 32) | (lower >> 32);
+    upper >>= 32;
+  }
+};
+
 class bigint {
  private:
   // A bigint is stored as an array of bigits (big digits), with bigit at index
@@ -517,7 +538,32 @@ class bigint {
   void operator=(const bigint&) = delete;
 
   void square() {
-    // TODO
+    basic_memory_buffer<bigit> n(std::move(bigits_));
+    int num_bigits = static_cast<int>(bigits_.size());
+    int num_result_bigints = 2 * num_bigits;
+    bigits_.resize(num_result_bigints);
+    using accumulator_t = conditional_t<FMT_USE_INT128, uint128_t, accumulator>;
+    auto sum = accumulator_t();
+    for (int bigit_index = 0; bigit_index < num_bigits; ++bigit_index) {
+      // Compute bigit at position bigit_index of the result by adding
+      // cross-product terms n[i] * n[j] such that i + j == bigit_index.
+      for (int i = 0, j = bigit_index; j >= 0; ++i, --j) {
+        // Most terms are multiplied twice which can be optimized in the future.
+        sum += static_cast<double_bigit>(n[i]) * n[j];
+      }
+      bigits_[bigit_index] = static_cast<bigit>(sum);
+      sum >>= bits<bigit>::value;  // Compute the carry.
+    }
+    // Do the same for the top half.
+    for (int bigit_index = num_bigits; bigit_index < num_result_bigints;
+         ++bigit_index) {
+      for (int j = num_bigits - 1, i = bigit_index - j; i < num_bigits;
+           ++i, --j) {
+        sum += static_cast<double_bigit>(n[i]) * n[j];
+      }
+      bigits_[bigit_index] = static_cast<bigit>(sum);
+      sum >>= bits<bigit>::value;
+    }
   }
 
   bigint& operator<<=(int shift) {
@@ -771,8 +817,12 @@ template <int GRISU_VERSION> struct grisu_shortest_handler {
 };
 
 FMT_FUNC void fallback_format(const fp& value, int exp10) {
-  (void)value; // TODO
-  (void)exp10;
+  bigint big_value(value.f), pow10;
+  if (value.e >= 0) {
+    big_value <<= value.e;
+    pow10.assign_pow10(exp10);
+  }
+  // TODO
 }
 
 template <typename Double,
@@ -869,12 +919,8 @@ char* sprintf_format(Double value, internal::buffer<char>& buf,
     type = 'f';
   else if (type == 0 || type == 'n')
     type = 'g';
-#if FMT_MSC_VER
-  if (type == 'F') {
-    // MSVC's printf doesn't support 'F'.
-    type = 'f';
-  }
-#endif
+  if (FMT_MSC_VER && type == 'F')
+    type = 'f';  // // MSVC's printf doesn't support 'F'.
   *format_ptr++ = type;
   *format_ptr = '\0';
 
diff --git a/test/format-impl-test.cc b/test/format-impl-test.cc
index 94db8d97..be81238e 100644
--- a/test/format-impl-test.cc
+++ b/test/format-impl-test.cc
@@ -59,6 +59,38 @@ TEST(BigIntTest, Multiply) {
   EXPECT_EQ("fffffffe00000001", fmt::format("{}", bigmax));
 }
 
+TEST(BigIntTest, Accumulator) {
+  fmt::internal::accumulator acc;
+  EXPECT_EQ(acc.lower, 0);
+  EXPECT_EQ(acc.upper, 0);
+  acc.upper = 12;
+  acc.lower = 34;
+  EXPECT_EQ(static_cast<uint32_t>(acc), 34);
+  acc += 56;
+  EXPECT_EQ(acc.lower, 90);
+  acc += fmt::internal::max_value<uint64_t>();
+  EXPECT_EQ(acc.upper, 13);
+  EXPECT_EQ(acc.lower, 89);
+  acc >>= 32;
+  EXPECT_EQ(acc.upper, 0);
+  EXPECT_EQ(acc.lower, 13 * 0x100000000);
+}
+
+TEST(BigIntTest, Square) {
+  bigint n0(0);
+  n0.square();
+  EXPECT_EQ("0", fmt::format("{}", n0));
+  bigint n1(0x100);
+  n1.square();
+  EXPECT_EQ("10000", fmt::format("{}", n1));
+  bigint n2(0xfffffffff);
+  n2.square();
+  EXPECT_EQ("ffffffffe000000001", fmt::format("{}", n2));
+  bigint n3(max_value<uint64_t>());
+  n3.square();
+  EXPECT_EQ("fffffffffffffffe0000000000000001", fmt::format("{}", n3));
+}
+
 template <bool is_iec559> void test_construct_from_double() {
   fmt::print("warning: double is not IEC559, skipping FP tests\n");
 }