[flang] Tuning up binary->decimal conversion

[klausler]

Summary

Use short division of big-radix values by powers of two when
converting values with negative unbiased exponents rather than
multiplication by smaller powers of five; this reduces the overall
outer iteration count. This change is a win across the entire range
of inputs.

flang/include/flang/Decimal/decimal.h

 
 /*
  * When allocating decimal conversion output buffers, use the maximum
  * number of significant decimal digits in the representation of the
  * least nonzero value, and add this extra space for a sign, a NUL, and
  * some extra due to the library working internally in base 10**16
  * and computing its output size in multiples of 16.
  */
-#define EXTRA_DECIMAL_CONVERSION_SPACE (1 + 1 + 16 - 1)
+#define EXTRA_DECIMAL_CONVERSION_SPACE (1 + 1 + 2 * 16 - 1)
 
 #ifdef __cplusplus
 template <int PREC>
 ConversionToDecimalResult ConvertToDecimal(char *, size_t,
     DecimalConversionFlags, int digits, enum FortranRounding rounding,
     BinaryFloatingPointNumber<PREC> x);
 
 extern template ConversionToDecimalResult ConvertToDecimal<8>(char *, size_t,

flang/lib/Decimal/big-radix-floating-point.h

for (int j{digits_ - 1}; j >= 0; --j) {
       Digit q{common::DivideUnsignedBy<Digit, DIVISOR>(digit_[j])};
       Digit nrem{digit_[j] - DIVISOR * q};
       digit_[j] = q + (radix / DIVISOR) * remainder;
       remainder = nrem;
     }
     return remainder;
   }
 
-  int DivideByPowerOfTwo(int twoPow) { // twoPow <= LOG10RADIX
-    int remainder{0};
+  void DivideByPowerOfTwo(int twoPow) { // twoPow <= log10Radix
+    Digit remainder{0};
+    auto mask{(Digit{1} << twoPow) - 1};
+    auto coeff{radix >> twoPow};
     for (int j{digits_ - 1}; j >= 0; --j) {
-      Digit q{digit_[j] >> twoPow};
-      int nrem = digit_[j] - (q << twoPow);
-      digit_[j] = q + (radix >> twoPow) * remainder;
+      auto nrem{digit_[j] & mask};
+      digit_[j] = (digit_[j] >> twoPow) + coeff * remainder;
       remainder = nrem;
     }
-    return remainder;
+  }
+
+  // Returns true on overflow
+  bool DivideByPowerOfTwoInPlace(int twoPow) {
+    if (digits_ > 0) {
+      while (twoPow > 0) {
+        int chunk{twoPow > log10Radix ? log10Radix : twoPow};
+        if ((digit_[0] & ((Digit{1} << chunk) - 1)) == 0) {
+          DivideByPowerOfTwo(chunk);
+          twoPow -= chunk;
+          continue;
+        }
+        twoPow -= chunk;
+        if (digit_[digits_ - 1] >> chunk != 0) {
+          if (digits_ == digitLimit_) {
+            return true; // overflow
+          }
+          digit_[digits_++] = 0;
+        }
+        auto remainder{digit_[digits_ - 1]};
+        exponent_ -= log10Radix;
+        auto coeff{radix >> chunk}; // precise; radix is (5*2)**log10Radix
+        auto mask{(Digit{1} << chunk) - 1};
+        for (int j{digits_ - 1}; j >= 1; --j) {
+          digit_[j] = (digit_[j - 1] >> chunk) + coeff * remainder;
+          remainder = digit_[j - 1] & mask;
+        }
+        digit_[0] = coeff * remainder;
+      }
+    }
+    return false; // no overflow
   }
 
   int AddCarry(int position = 0, int carry = 1) {
     for (; position < digits_; ++position) {
       Digit v{digit_[position] + carry};
       if (v < radix) {
         digit_[position] = v;
         return 0;

flang/lib/Decimal/binary-to-decimal.cpp

for (; twoPow >= 3; twoPow -= 3) {
     // D * 10.**E * 2.**twoPow -> (D*(2**3)) * 10.**E * 2.**(twoPow-3)
     overflow |= MultiplyBy<8>();
   }
   for (; twoPow > 0; --twoPow) {
     // D * 10.**E * 2.**twoPow -> (2*D) * 10.**E * 2.**(twoPow-1)
     overflow |= MultiplyBy<2>();
   }
 
-  while (twoPow < 0) {
-    int shift{common::TrailingZeroBitCount(digit_[0])};
-    if (shift == 0) {
-      break;
-    }
-    if (shift > log10Radix) {
-      shift = log10Radix;
-    }
-    if (shift > -twoPow) {
-      shift = -twoPow;
-    }
-    // (D*(2**S)) * 10.**E * 2.**twoPow -> D * 10.**E * 2.**(twoPow+S)
-    DivideByPowerOfTwo(shift);
-    twoPow += shift;
-  }
-
-  for (; twoPow <= -4; twoPow += 4) {
-    // D * 10.**E * 2.**twoPow -> 625D * 10.**(E-4) * 2.**(twoPow+4)
-    overflow |= MultiplyBy<(5 * 5 * 5 * 5)>();
-    exponent_ -= 4;
-  }
-  if (twoPow <= -2) {
-    // D * 10.**E * 2.**twoPow -> 25D * 10.**(E-2) * 2.**(twoPow+2)
-    overflow |= MultiplyBy<5 * 5>();
-    twoPow += 2;
-    exponent_ -= 2;
-  }
-  for (; twoPow < 0; ++twoPow) {
-    // D * 10.**E * 2.**twoPow -> 5D * 10.**(E-1) * 2.**(twoPow+1)
-    overflow |= MultiplyBy<5>();
-    --exponent_;
-  }
-
+  overflow |= DivideByPowerOfTwoInPlace(-twoPow);
   assert(overflow == 0);
-
-  // twoPow == 0, the decimal encoding is complete.
   Normalize();
 }
 
 template <int PREC, int LOG10RADIX>
 ConversionToDecimalResult
 BigRadixFloatingPointNumber<PREC, LOG10RADIX>::ConvertToDecimal(char *buffer,
     std::size_t n, enum DecimalConversionFlags flags, int maxDigits) const {
   if (n < static_cast<std::size_t>(3 + digits_ * LOG10RADIX)) {
if (q[0] != '0' || p > start) {
       *p++ = q[1];
     }
   }
   for (int j{digits_ - 1}; j-- > 0;) {
     Digit dig{digit_[j]};
     for (int k{0}; k < log10Radix; k += 2) {
       Digit d{common::DivideUnsignedBy<Digit, hundredth>(dig)};
       dig = 100 * (dig - d * hundredth);
-      const char *q = lut + 2 * d;
+      const char *q{lut + 2 * d};
       *p++ = q[0];
       *p++ = q[1];
     }
   }
   // Adjust exponent so the effective decimal point is to
   // the left of the first digit.
   int expo = exponent_ + p - start;
   // Trim trailing zeroes.