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 <cfloat>
#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()};
  if (x.IsZero()) {
    isNegative_ = negative;
    return;
  }
  if (negative) {
    x.Negate();
  }
  int twoPow{x.UnbiasedExponent()};
  twoPow -= x.bits - 1;
  if (!x.isImplicitMSB) {
    ++twoPow;
  }
  int lshift{x.exponentBits};
  if (twoPow <= -lshift) {
    twoPow += lshift;
    lshift = 0;
  } else if (twoPow < 0) {
    lshift += twoPow;
    twoPow = 0;
  }
  auto word{x.Fraction()};
  word <<= lshift;
  SetTo(word);
  isNegative_ = negative;

  // The significand is now encoded in *this as an integer (D) and
  // decimal exponent (E):  x = D * 10.**E * 2.**twoPow
  // twoPow can be positive or negative.
  // The goal now is to get twoPow up or down to zero, leaving us with
  // only decimal digits and decimal exponent.  This is done by
  // 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)
    overflow |= MultiplyBy<512>();
  }
  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>();
  }

  overflow |= DivideByPowerOfTwoInPlace(-twoPow);
  assert(overflow == 0);
  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)) {
    return {nullptr, 0, 0, Overflow};
  }
  char *start{buffer};
  if (isNegative_) {
    *start++ = '-';
  } else if (flags & AlwaysSign) {
    *start++ = '+';
  }
  if (IsZero()) {
    *start++ = '0';
    *start = '\0';
    return {buffer, static_cast<std::size_t>(start - buffer), 0, Exact};
  }
  char *p{start};
  static_assert((LOG10RADIX % 2) == 0, "radix not a power of 100");
  static const char lut[] = "0001020304050607080910111213141516171819"
                            "2021222324252627282930313233343536373839"
                            "4041424344454647484950515253545556575859"
                            "6061626364656667686970717273747576777879"
                            "8081828384858687888990919293949596979899";
  // Treat the MSD specially: don't emit leading zeroes.
  Digit dig{digit_[digits_ - 1]};
  char stack[LOG10RADIX], *sp{stack};
  for (int k{0}; k < log10Radix; k += 2) {
    Digit newDig{dig / 100};
    auto d{static_cast<std::uint32_t>(dig) -
        std::uint32_t{100} * static_cast<std::uint32_t>(newDig)};
    dig = newDig;
    const char *q{lut + d + d};
    *sp++ = q[1];
    *sp++ = q[0];
  }
  while (sp > stack && sp[-1] == '0') {
    --sp;
  }
  while (sp > stack) {
    *p++ = *--sp;
  }
  for (int j{digits_ - 1}; j-- > 0;) {
    Digit dig{digit_[j]};
    char *reverse{p += log10Radix};
    for (int k{0}; k < log10Radix; k += 2) {
      Digit newDig{dig / 100};
      auto d{static_cast<std::uint32_t>(dig) -
          std::uint32_t{100} * static_cast<std::uint32_t>(newDig)};
      dig = newDig;
      const char *q{lut + d + d};
      *--reverse = q[1];
      *--reverse = q[0];
    }
  }
  // Adjust exponent so the effective decimal point is to
  // the left of the first digit.
  int expo = exponent_ + p - start;
  // Trim trailing zeroes.
  while (p[-1] == '0') {
    --p;
  }
  char *end{start + maxDigits};
  if (maxDigits == 0) {
    p = end;
  }
  if (p <= end) {
    *p = '\0';
    return {buffer, static_cast<std::size_t>(p - buffer), expo, Exact};
  } else {
    // Apply a digit limit, possibly with rounding.
    bool incr{false};
    switch (rounding_) {
    case RoundNearest:
      incr = *end > '5' ||
          (*end == '5' && (p > end + 1 || ((end[-1] - '0') & 1) != 0));
      break;
    case RoundUp:
      incr = !isNegative_;
      break;
    case RoundDown:
      incr = isNegative_;
      break;
    case RoundToZero:
      break;
    case RoundCompatible:
      incr = *end >= '5';
      break;
    }
    p = end;
    if (incr) {
      while (p > start && p[-1] == '9') {
        --p;
      }
      if (p == start) {
        *p++ = '1';
        ++expo;
      } else {
        ++p[-1];
      }
    }

    *p = '\0';
    return {buffer, static_cast<std::size_t>(p - buffer), expo, Inexact};
  }
}

template <int PREC, int LOG10RADIX>
bool BigRadixFloatingPointNumber<PREC, LOG10RADIX>::Mean(
    const BigRadixFloatingPointNumber &that) {
  while (digits_ < that.digits_) {
    digit_[digits_++] = 0;
  }
  int carry{0};
  for (int j{0}; j < that.digits_; ++j) {
    Digit v{digit_[j] + that.digit_[j] + carry};
    if (v >= radix) {
      digit_[j] = v - radix;
      carry = 1;
    } else {
      digit_[j] = v;
      carry = 0;
    }
  }
  if (carry != 0) {
    AddCarry(that.digits_, carry);
  }
  return DivideBy<2>() != 0;
}

template <int PREC, int LOG10RADIX>
void BigRadixFloatingPointNumber<PREC, LOG10RADIX>::Minimize(
    BigRadixFloatingPointNumber &&less, BigRadixFloatingPointNumber &&more) {
  int leastExponent{exponent_};
  if (less.exponent_ < leastExponent) {
    leastExponent = less.exponent_;
  }
  if (more.exponent_ < leastExponent) {
    leastExponent = more.exponent_;
  }
  while (exponent_ > leastExponent) {
    --exponent_;
    MultiplyBy<10>();
  }
  while (less.exponent_ > leastExponent) {
    --less.exponent_;
    less.MultiplyBy<10>();
  }
  while (more.exponent_ > leastExponent) {
    --more.exponent_;
    more.MultiplyBy<10>();
  }
  if (less.Mean(*this)) {
    less.AddCarry(); // round up
  }
  if (!more.Mean(*this)) {
    more.Decrement(); // round down
  }
  while (less.digits_ < more.digits_) {
    less.digit_[less.digits_++] = 0;
  }
  while (more.digits_ < less.digits_) {
    more.digit_[more.digits_++] = 0;
  }
  int digits{more.digits_};
  int same{0};
  while (same < digits &&
      less.digit_[digits - 1 - same] == more.digit_[digits - 1 - same]) {
    ++same;
  }
  if (same == digits) {
    return;
  }
  digits_ = same + 1;
  int offset{digits - digits_};
  exponent_ += offset * log10Radix;
  for (int j{0}; j < digits_; ++j) {
    digit_[j] = more.digit_[j + offset];
  }
  Digit least{less.digit_[offset]};
  Digit my{digit_[0]};
  while (true) {
    Digit q{my / 10u};
    Digit r{my - 10 * q};
    Digit lq{least / 10u};
    Digit lr{least - 10 * lq};
    if (r != 0 && lq == q) {
      Digit sub{(r - lr) >> 1};
      digit_[0] -= sub;
      break;
    } else {
      least = lq;
      my = q;
      DivideBy<10>();
      ++exponent_;
    }
  }
  Normalize();
}

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()) {
      return {"-Inf", 4, 0, Exact};
    } else if (flags & AlwaysSign) {
      return {"+Inf", 4, 0, Exact};
    } else {
      return {"Inf", 3, 0, Exact};
    }
  } else {
    using Big = BigRadixFloatingPointNumber<PREC>;
    Big number{x, rounding};
    if ((flags & Minimize) && !x.IsZero()) {
      // To emit the fewest decimal digits necessary to represent the value
      // in such a way that decimal-to-binary conversion to the same format
      // 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.Previous();
      Binary more{x};
      if (!x.IsMaximalFiniteMagnitude()) {
        more.Next();
      }
      number.Minimize(Big{less, rounding}, Big{more, rounding});
    }
    return number.ConvertToDecimal(buffer, size, flags, digits);
  }
}

template ConversionToDecimalResult ConvertToDecimal<8>(char *, std::size_t,
    enum DecimalConversionFlags, int, enum FortranRounding,
    BinaryFloatingPointNumber<8>);
template ConversionToDecimalResult ConvertToDecimal<11>(char *, std::size_t,
    enum DecimalConversionFlags, int, enum FortranRounding,
    BinaryFloatingPointNumber<11>);
template ConversionToDecimalResult ConvertToDecimal<24>(char *, std::size_t,
    enum DecimalConversionFlags, int, enum FortranRounding,
    BinaryFloatingPointNumber<24>);
template ConversionToDecimalResult ConvertToDecimal<53>(char *, std::size_t,
    enum DecimalConversionFlags, int, enum FortranRounding,
    BinaryFloatingPointNumber<53>);
template ConversionToDecimalResult ConvertToDecimal<64>(char *, std::size_t,
    enum DecimalConversionFlags, int, enum FortranRounding,
    BinaryFloatingPointNumber<64>);
template ConversionToDecimalResult ConvertToDecimal<113>(char *, std::size_t,
    enum DecimalConversionFlags, int, enum FortranRounding,
    BinaryFloatingPointNumber<113>);

extern "C" {
RT_EXT_API_GROUP_BEGIN

ConversionToDecimalResult ConvertFloatToDecimal(char *buffer, std::size_t size,
    enum DecimalConversionFlags flags, int digits,
    enum FortranRounding rounding, float x) {
  return Fortran::decimal::ConvertToDecimal(buffer, size, flags, digits,
      rounding, Fortran::decimal::BinaryFloatingPointNumber<24>(x));
}

ConversionToDecimalResult ConvertDoubleToDecimal(char *buffer, std::size_t size,
    enum DecimalConversionFlags flags, int digits,
    enum FortranRounding rounding, double x) {
  return Fortran::decimal::ConvertToDecimal(buffer, size, flags, digits,
      rounding, Fortran::decimal::BinaryFloatingPointNumber<53>(x));
}

#if LDBL_MANT_DIG == 64
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));
}
#elif LDBL_MANT_DIG == 113
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<113>(x));
}
#endif

RT_EXT_API_GROUP_END
} // extern "C"

template <int PREC, int LOG10RADIX>
template <typename STREAM>
STREAM &BigRadixFloatingPointNumber<PREC, LOG10RADIX>::Dump(STREAM &o) const {
  if (isNegative_) {
    o << '-';
  }
  o << "10**(" << exponent_ << ") * ...  (rounding "
    << static_cast<int>(rounding_) << ")\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