[flang] Fix x87 binary->decimal

[klausler]

Summary

Fix decimal formatting of 80-bit x87 values; the calculation ofnearest neighbor values failed to account for the explicitmost significant bit in that format.

Replace MultiplyByRounded with MultiplyBy in binary->decimal conversions,
since rounding won't happen and the name was misleading; then remove
dead code, and migrate LoseLeastSignificantDigit() from one source file
to another where it's still needed.

flang/include/flang/Common/uint128.h

 //===-- include/flang/Common/uint128.h --------------------------*- C++ -*-===//
 //
 // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
 // See https://llvm.org/LICENSE.txt for license information.
 // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
 //
 //===----------------------------------------------------------------------===//
 
 // Portable 128-bit unsigned integer arithmetic for use in impoverished
 // C++ implementations lacking __uint128_t.
 
 #ifndef FORTRAN_COMMON_UINT128_H_
 #define FORTRAN_COMMON_UINT128_H_
 
+// Define AVOID_NATIVE_UINT128_T to force the use of UnsignedInt128 below
+// instead of the C++ compiler's native 128-bit unsigned integer type, if
+// it has one.
 #ifndef AVOID_NATIVE_UINT128_T
-#define AVOID_NATIVE_UINT128_T 1 // always use this code for now for testing
+#define AVOID_NATIVE_UINT128_T 0
 #endif
 
 #include "leading-zero-bit-count.h"
 #include <cstdint>
 #include <type_traits>
 
 namespace Fortran::common {
 

flang/include/flang/Decimal/binary-floating-point.h

 #include <cinttypes>
 #include <climits>
 #include <cstring>
 #include <type_traits>
 
 namespace Fortran::decimal {
 
 template <int BINARY_PRECISION>
-struct BinaryFloatingPointNumber
-    : public common::RealDetails<BINARY_PRECISION> {
-
+class BinaryFloatingPointNumber : public common::RealDetails<BINARY_PRECISION> {
+public:
   using Details = common::RealDetails<BINARY_PRECISION>;
   using Details::bits;
   using Details::decimalPrecision;
   using Details::decimalRange;
   using Details::exponentBias;
   using Details::exponentBits;
   using Details::isImplicitMSB;
   using Details::maxDecimalConversionDigits;
constexpr BinaryFloatingPointNumber(
       const BinaryFloatingPointNumber &that) = default;
   constexpr BinaryFloatingPointNumber(
       BinaryFloatingPointNumber &&that) = default;
   constexpr BinaryFloatingPointNumber &operator=(
       const BinaryFloatingPointNumber &that) = default;
   constexpr BinaryFloatingPointNumber &operator=(
       BinaryFloatingPointNumber &&that) = default;
 
+  RawType raw() const { return raw_; }
+
   template <typename A> explicit constexpr BinaryFloatingPointNumber(A x) {
-    static_assert(sizeof raw <= sizeof x);
-    std::memcpy(reinterpret_cast<void *>(&raw),
-        reinterpret_cast<const void *>(&x), sizeof raw);
+    static_assert(sizeof raw_ <= sizeof x);
+    std::memcpy(reinterpret_cast<void *>(&raw_),
+        reinterpret_cast<const void *>(&x), sizeof raw_);
   }
 
   constexpr int BiasedExponent() const {
     return static_cast<int>(
-        (raw >> significandBits) & ((1 << exponentBits) - 1));
+        (raw_ >> significandBits) & ((1 << exponentBits) - 1));
   }
   constexpr int UnbiasedExponent() const {
     int biased{BiasedExponent()};
     return biased - exponentBias + (biased == 0);
   }
-  constexpr RawType Significand() const { return raw & significandMask; }
+  constexpr RawType Significand() const { return raw_ & significandMask; }
   constexpr RawType Fraction() const {
     RawType sig{Significand()};
     if (isImplicitMSB && BiasedExponent() > 0) {
       sig |= RawType{1} << significandBits;
     }
     return sig;
   }
 
   constexpr bool IsZero() const {
-    return (raw & ((RawType{1} << (bits - 1)) - 1)) == 0;
+    return (raw_ & ((RawType{1} << (bits - 1)) - 1)) == 0;
   }
   constexpr bool IsNaN() const {
     return BiasedExponent() == maxExponent && Significand() != 0;
   }
   constexpr bool IsInfinite() const {
     return BiasedExponent() == maxExponent && Significand() == 0;
   }
   constexpr bool IsMaximalFiniteMagnitude() const {
     return BiasedExponent() == maxExponent - 1 &&
         Significand() == significandMask;
   }
-  constexpr bool IsNegative() const { return ((raw >> (bits - 1)) & 1) != 0; }
+  constexpr bool IsNegative() const { return ((raw_ >> (bits - 1)) & 1) != 0; }
+
+  constexpr void Negate() { raw_ ^= RawType{1} << (bits - 1); }
+
+  // For calculating the nearest neighbors of a floating-point value
+  constexpr void Previous() {
+    RemoveExplicitMSB();
+    --raw_;
+    InsertExplicitMSB();
+  }
+  constexpr void Next() {
+    RemoveExplicitMSB();
+    ++raw_;
+    InsertExplicitMSB();
+  }
 
-  constexpr void Negate() { raw ^= RawType{1} << (bits - 1); }
+private:
+  constexpr void RemoveExplicitMSB() {
+    if constexpr (!isImplicitMSB) {
+      raw_ = (raw_ & (significandMask >> 1)) | ((raw_ & ~significandMask) >> 1);
+    }
+  }
+  constexpr void InsertExplicitMSB() {
+    if constexpr (!isImplicitMSB) {
+      constexpr RawType mask{significandMask >> 1};
+      raw_ = (raw_ & mask) | ((raw_ & ~mask) << 1);
+      if (BiasedExponent() > 0) {
+        raw_ |= RawType{1} << (significandBits - 1);
+      }
+    }
+  }
 
-  RawType raw{0};
+  RawType raw_{0};
 };
 } // namespace Fortran::decimal
 #endif

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

 // a general-purpose facility.
 
 #include "flang/Common/bit-population-count.h"
 #include "flang/Common/leading-zero-bit-count.h"
 #include "flang/Common/uint128.h"
 #include "flang/Common/unsigned-const-division.h"
 #include "flang/Decimal/binary-floating-point.h"
 #include "flang/Decimal/decimal.h"
+#include "llvm/Support/raw_ostream.h"
 #include <cinttypes>
 #include <limits>
 #include <type_traits>
 
 namespace Fortran::decimal {
 
 static constexpr std::uint64_t TenToThe(int power) {
   return power <= 0 ? 1 : 10 * TenToThe(power - 1);
public:
   // cases, and will be the mean value when that is so (e.g., if
   // last decimal digit values 6-8 would all work, it'll be a 7).
   // This minimization necessarily assumes that the value will be
   // emitted and read back into the same (or less precise) format
   // with default rounding to the nearest value.
   void Minimize(
       BigRadixFloatingPointNumber &&less, BigRadixFloatingPointNumber &&more);
 
+  llvm::raw_ostream &Dump(llvm::raw_ostream &) const;
+
 private:
   BigRadixFloatingPointNumber(const BigRadixFloatingPointNumber &that)
       : digits_{that.digits_}, exponent_{that.exponent_},
         isNegative_{that.isNegative_}, rounding_{that.rounding_} {
     for (int j{0}; j < digits_; ++j) {
       digit_[j] = that.digit_[j];
     }
   }
if (int carry{MultiplyByHelper<N>(0)}) {
       } else {
         return carry;
       }
     } else {
       return 0;
     }
   }
 
-  template <int N> void MultiplyByRounded() {
-    if (int carry{MultiplyBy<N>()}) {
-      LoseLeastSignificantDigit();
-      digit_[digits_ - 1] += carry;
-      exponent_ += log10Radix;
-    }
-  }
-
   void LoseLeastSignificantDigit(); // with rounding
 
   void PushCarry(int carry) {
     if (digits_ == maxDigits && RemoveLeastOrderZeroDigits() == 0) {
       LoseLeastSignificantDigit();
       digit_[digits_ - 1] += carry;
     } else {
       digit_[digits_++] = carry;

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

 //===-- lib/Decimal/binary-to-decimal.cpp ---------------------------------===//
 //
 // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
 // See https://llvm.org/LICENSE.txt for license information.
 // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
 //
 //===----------------------------------------------------------------------===//
 
 #include "big-radix-floating-point.h"
 #include "flang/Decimal/decimal.h"
+#include <cassert>
+#include <string>
 
 namespace Fortran::decimal {
 
 template <int PREC, int LOG10RADIX>
 BigRadixFloatingPointNumber<PREC, LOG10RADIX>::BigRadixFloatingPointNumber(
     BinaryFloatingPointNumber<PREC> x, enum FortranRounding rounding)
     : rounding_{rounding} {
   bool negative{x.IsNegative()};
BigRadixFloatingPointNumber<PREC, LOG10RADIX>::BigRadixFloatingPointNumber(
   // fast multiplications and divisions of D by 2 and 5.
 
   // (5*D) * 10.**E * 2.**twoPow -> D * 10.**(E+1) * 2.**(twoPow-1)
   for (; twoPow > 0 && IsDivisibleBy<5>(); --twoPow) {
     DivideBy<5>();
     ++exponent_;
   }
 
+  int overflow{0};
   for (; twoPow >= 9; twoPow -= 9) {
     // D * 10.**E * 2.**twoPow -> (D*(2**9)) * 10.**E * 2.**(twoPow-9)
-    MultiplyByRounded<512>();
+    overflow |= MultiplyBy<512>();
   }
   for (; twoPow >= 3; twoPow -= 3) {
     // D * 10.**E * 2.**twoPow -> (D*(2**3)) * 10.**E * 2.**(twoPow-3)
-    MultiplyByRounded<8>();
+    overflow |= MultiplyBy<8>();
   }
   for (; twoPow > 0; --twoPow) {
     // D * 10.**E * 2.**twoPow -> (2*D) * 10.**E * 2.**(twoPow-1)
-    MultiplyByRounded<2>();
+    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)
-    MultiplyByRounded<(5 * 5 * 5 * 5)>();
+    overflow |= MultiplyBy<(5 * 5 * 5 * 5)>();
     exponent_ -= 4;
   }
   if (twoPow <= -2) {
     // D * 10.**E * 2.**twoPow -> 25D * 10.**(E-2) * 2.**(twoPow+2)
-    MultiplyByRounded<25>();
+    overflow |= MultiplyBy<5 * 5>();
     twoPow += 2;
     exponent_ -= 2;
   }
   for (; twoPow < 0; ++twoPow) {
     // D * 10.**E * 2.**twoPow -> 5D * 10.**(E-1) * 2.**(twoPow+1)
-    MultiplyByRounded<5>();
+    overflow |= MultiplyBy<5>();
     --exponent_;
   }
 
+  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 (r != 0 && lq == q) {
       my = q;
       DivideBy<10>();
       ++exponent_;
     }
   }
   Normalize();
 }
 
-template <int PREC, int LOG10RADIX>
-void BigRadixFloatingPointNumber<PREC,
-    LOG10RADIX>::LoseLeastSignificantDigit() {
-  Digit LSD{digit_[0]};
-  for (int j{0}; j < digits_ - 1; ++j) {
-    digit_[j] = digit_[j + 1];
-  }
-  digit_[digits_ - 1] = 0;
-  bool incr{false};
-  switch (rounding_) {
-  case RoundNearest:
-  case RoundDefault:
-    incr = LSD > radix / 2 || (LSD == radix / 2 && digit_[0] % 2 != 0);
-    break;
-  case RoundUp:
-    incr = LSD > 0 && !isNegative_;
-    break;
-  case RoundDown:
-    incr = LSD > 0 && isNegative_;
-    break;
-  case RoundToZero:
-    break;
-  case RoundCompatible:
-    incr = LSD >= radix / 2;
-    break;
-  }
-  for (int j{0}; (digit_[j] += incr) == radix; ++j) {
-    digit_[j] = 0;
-  }
-}
-
 template <int PREC>
 ConversionToDecimalResult ConvertToDecimal(char *buffer, std::size_t size,
     enum DecimalConversionFlags flags, int digits,
     enum FortranRounding rounding, BinaryFloatingPointNumber<PREC> x) {
   if (x.IsNaN()) {
     return {"NaN", 3, 0, Invalid};
   } else if (x.IsInfinite()) {
     if (x.IsNegative()) {
if ((flags & Minimize) && !x.IsZero()) {
       // with a fixed assumption about rounding will return the same binary
       // value, we also perform binary-to-decimal conversion on the two
       // binary values immediately adjacent to this one, use them to identify
       // the bounds of the range of decimal values that will map back to the
       // original binary value, and find a (not necessary unique) shortest
       // decimal sequence in that range.
       using Binary = typename Big::Real;
       Binary less{x};
-      --less.raw;
+      less.Previous();
       Binary more{x};
       if (!x.IsMaximalFiniteMagnitude()) {
-        ++more.raw;
+        more.Next();
       }
       number.Minimize(Big{less, rounding}, Big{more, rounding});
+    } else {
     }
     return number.ConvertToDecimal(buffer, size, flags, digits);
   }
 }
 
 template ConversionToDecimalResult ConvertToDecimal<8>(char *, std::size_t,
     enum DecimalConversionFlags, int, enum FortranRounding,
     BinaryFloatingPointNumber<8>);

 ConversionToDecimalResult ConvertLongDoubleToDecimal(char *buffer,
     std::size_t size, enum DecimalConversionFlags flags, int digits,
     enum FortranRounding rounding, long double x) {
   return Fortran::decimal::ConvertToDecimal(buffer, size, flags, digits,
       rounding, Fortran::decimal::BinaryFloatingPointNumber<64>(x));
 }
 #endif
 }
+
+template <int PREC, int LOG10RADIX>
+llvm::raw_ostream &BigRadixFloatingPointNumber<PREC, LOG10RADIX>::Dump(
+    llvm::raw_ostream &o) const {
+  if (isNegative_) {
+    o << '-';
+  }
+  o << "10**(" << exponent_ << ") * ...\n";
+  for (int j{digits_}; --j >= 0;) {
+    std::string str{std::to_string(digit_[j])};
+    o << std::string(20 - str.size(), ' ') << str << " [" << j << ']';
+    if (j + 1 == digitLimit_) {
+      o << " (limit)";
+    }
+    o << '\n';
+  }
+  return o;
+}
 } // namespace Fortran::decimal

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

case 'Q': {
     }
   } break;
   default:
     break;
   }
   return true;
 }
 
+template <int PREC, int LOG10RADIX>
+void BigRadixFloatingPointNumber<PREC,
+    LOG10RADIX>::LoseLeastSignificantDigit() {
+  Digit LSD{digit_[0]};
+  for (int j{0}; j < digits_ - 1; ++j) {
+    digit_[j] = digit_[j + 1];
+  }
+  digit_[digits_ - 1] = 0;
+  bool incr{false};
+  switch (rounding_) {
+  case RoundNearest:
+  case RoundDefault:
+    incr = LSD > radix / 2 || (LSD == radix / 2 && digit_[0] % 2 != 0);
+    break;
+  case RoundUp:
+    incr = LSD > 0 && !isNegative_;
+    break;
+  case RoundDown:
+    incr = LSD > 0 && isNegative_;
+    break;
+  case RoundToZero:
+    break;
+  case RoundCompatible:
+    incr = LSD >= radix / 2;
+    break;
+  }
+  for (int j{0}; (digit_[j] += incr) == radix; ++j) {
+    digit_[j] = 0;
+  }
+}
+
 // This local utility class represents an unrounded nonnegative
 // binary floating-point value with an unbiased (i.e., signed)
 // binary exponent, an integer value (not a fraction) with an implied
 // binary point to its *right*, and some guard bits for rounding.
 template <int PREC> class IntermediateFloat {
 public:
   static constexpr int precision{PREC};
   using IntType = common::HostUnsignedIntType<precision>;

flang/runtime/edit-output.cpp

case 'E':
     } else {
       return EditEorDOutput(edit);
     }
   case 'F':
     return EditFOutput(edit);
   case 'B':
   case 'O':
   case 'Z':
-    return EditIntegerOutput(
-        io_, edit, decimal::BinaryFloatingPointNumber<binaryPrecision>{x_}.raw);
+    return EditIntegerOutput(io_, edit,
+        decimal::BinaryFloatingPointNumber<binaryPrecision>{x_}.raw());
   case 'G':
     return Edit(EditForGOutput(edit));
   default:
     if (edit.IsListDirected()) {
       return EditListDirectedOutput(edit);
     }
     io_.GetIoErrorHandler().SignalError(IostatErrorInFormat,
         "Data edit descriptor '%c' may not be used with a REAL data item",