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

   1 //===-- lib/Decimal/binary-to-decimal.cpp ---------------------------------===//
   2 //
   3 // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
   4 // See https://llvm.org/LICENSE.txt for license information.
   5 // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
   6 //
   7 //===----------------------------------------------------------------------===//
   8
   9 #include "big-radix-floating-point.h"
  10 #include "flang/Decimal/decimal.h"
  11 #include <cassert>
  12 #include <cfloat>
  13 #include <string>
  14
  15 namespace Fortran::decimal {
  16
  17 template <int PREC, int LOG10RADIX>
  18 BigRadixFloatingPointNumber<PREC, LOG10RADIX>::BigRadixFloatingPointNumber(
  19     BinaryFloatingPointNumber<PREC> x, enum FortranRounding rounding)
  20     : rounding_{rounding} {
  21   bool negative{x.IsNegative()};
  22   if (x.IsZero()) {
  23     isNegative_ = negative;
  24     return;
  25   }
  26   if (negative) {
  27     x.Negate();
  28   }
  29   int twoPow{x.UnbiasedExponent()};
  30   twoPow -= x.bits - 1;
  31   if (!x.isImplicitMSB) {
  32     ++twoPow;
  33   }
  34   int lshift{x.exponentBits};
  35   if (twoPow <= -lshift) {
  36     twoPow += lshift;
  37     lshift = 0;
  38   } else if (twoPow < 0) {
  39     lshift += twoPow;
  40     twoPow = 0;
  41   }
  42   auto word{x.Fraction()};
  43   word <<= lshift;
  44   SetTo(word);
  45   isNegative_ = negative;
  46
  47   // The significand is now encoded in *this as an integer (D) and
  48   // decimal exponent (E):  x = D * 10.**E * 2.**twoPow
  49   // twoPow can be positive or negative.
  50   // The goal now is to get twoPow up or down to zero, leaving us with
  51   // only decimal digits and decimal exponent.  This is done by
  52   // fast multiplications and divisions of D by 2 and 5.
  53
  54   // (5*D) * 10.**E * 2.**twoPow -> D * 10.**(E+1) * 2.**(twoPow-1)
  55   for (; twoPow > 0 && IsDivisibleBy<5>(); --twoPow) {
  56     DivideBy<5>();
  57     ++exponent_;
  58   }
  59
  60   int overflow{0};
  61   for (; twoPow >= 9; twoPow -= 9) {
  62     // D * 10.**E * 2.**twoPow -> (D*(2**9)) * 10.**E * 2.**(twoPow-9)
  63     overflow |= MultiplyBy<512>();
  64   }
  65   for (; twoPow >= 3; twoPow -= 3) {
  66     // D * 10.**E * 2.**twoPow -> (D*(2**3)) * 10.**E * 2.**(twoPow-3)
  67     overflow |= MultiplyBy<8>();
  68   }
  69   for (; twoPow > 0; --twoPow) {
  70     // D * 10.**E * 2.**twoPow -> (2*D) * 10.**E * 2.**(twoPow-1)
  71     overflow |= MultiplyBy<2>();
  72   }
  73
  74   overflow |= DivideByPowerOfTwoInPlace(-twoPow);
  75   assert(overflow == 0);
  76   Normalize();
  77 }
  78
  79 template <int PREC, int LOG10RADIX>
  80 ConversionToDecimalResult
  81 BigRadixFloatingPointNumber<PREC, LOG10RADIX>::ConvertToDecimal(char *buffer,
  82     std::size_t n, enum DecimalConversionFlags flags, int maxDigits) const {
  83   if (n < static_cast<std::size_t>(3 + digits_ * LOG10RADIX)) {
  84     return {nullptr, 0, 0, Overflow};
  85   }
  86   char *start{buffer};
  87   if (isNegative_) {
  88     *start++ = '-';
  89   } else if (flags & AlwaysSign) {
  90     *start++ = '+';
  91   }
  92   if (IsZero()) {
  93     *start++ = '0';
  94     *start = '\0';
  95     return {buffer, static_cast<std::size_t>(start - buffer), 0, Exact};
  96   }
  97   char *p{start};
  98   static_assert((LOG10RADIX % 2) == 0, "radix not a power of 100");
  99   static const char lut[] = "0001020304050607080910111213141516171819"
 100                             "2021222324252627282930313233343536373839"
 101                             "4041424344454647484950515253545556575859"
 102                             "6061626364656667686970717273747576777879"
 103                             "8081828384858687888990919293949596979899";
 104   // Treat the MSD specially: don't emit leading zeroes.
 105   Digit dig{digit_[digits_ - 1]};
 106   char stack[LOG10RADIX], *sp{stack};
 107   for (int k{0}; k < log10Radix; k += 2) {
 108     Digit newDig{dig / 100};
 109     auto d{static_cast<std::uint32_t>(dig) -
 110         std::uint32_t{100} * static_cast<std::uint32_t>(newDig)};
 111     dig = newDig;
 112     const char *q{lut + d + d};
 113     *sp++ = q[1];
 114     *sp++ = q[0];
 115   }
 116   while (sp > stack && sp[-1] == '0') {
 117     --sp;
 118   }
 119   while (sp > stack) {
 120     *p++ = *--sp;
 121   }
 122   for (int j{digits_ - 1}; j-- > 0;) {
 123     Digit dig{digit_[j]};
 124     char *reverse{p += log10Radix};
 125     for (int k{0}; k < log10Radix; k += 2) {
 126       Digit newDig{dig / 100};
 127       auto d{static_cast<std::uint32_t>(dig) -
 128           std::uint32_t{100} * static_cast<std::uint32_t>(newDig)};
 129       dig = newDig;
 130       const char *q{lut + d + d};
 131       *--reverse = q[1];
 132       *--reverse = q[0];
 133     }
 134   }
 135   // Adjust exponent so the effective decimal point is to
 136   // the left of the first digit.
 137   int expo = exponent_ + p - start;
 138   // Trim trailing zeroes.
 139   while (p[-1] == '0') {
 140     --p;
 141   }
 142   char *end{start + maxDigits};
 143   if (maxDigits == 0) {
 144     p = end;
 145   }
 146   if (p <= end) {
 147     *p = '\0';
 148     return {buffer, static_cast<std::size_t>(p - buffer), expo, Exact};
 149   } else {
 150     // Apply a digit limit, possibly with rounding.
 151     bool incr{false};
 152     switch (rounding_) {
 153     case RoundNearest:
 154       incr = *end > '5' ||
 155           (*end == '5' && (p > end + 1 || ((end[-1] - '0') & 1) != 0));
 156       break;
 157     case RoundUp:
 158       incr = !isNegative_;
 159       break;
 160     case RoundDown:
 161       incr = isNegative_;
 162       break;
 163     case RoundToZero:
 164       break;
 165     case RoundCompatible:
 166       incr = *end >= '5';
 167       break;
 168     }
 169     p = end;
 170     if (incr) {
 171       while (p > start && p[-1] == '9') {
 172         --p;
 173       }
 174       if (p == start) {
 175         *p++ = '1';
 176         ++expo;
 177       } else {
 178         ++p[-1];
 179       }
 180     }
 181
 182     *p = '\0';
 183     return {buffer, static_cast<std::size_t>(p - buffer), expo, Inexact};
 184   }
 185 }
 186
 187 template <int PREC, int LOG10RADIX>
 188 bool BigRadixFloatingPointNumber<PREC, LOG10RADIX>::Mean(
 189     const BigRadixFloatingPointNumber &that) {
 190   while (digits_ < that.digits_) {
 191     digit_[digits_++] = 0;
 192   }
 193   int carry{0};
 194   for (int j{0}; j < that.digits_; ++j) {
 195     Digit v{digit_[j] + that.digit_[j] + carry};
 196     if (v >= radix) {
 197       digit_[j] = v - radix;
 198       carry = 1;
 199     } else {
 200       digit_[j] = v;
 201       carry = 0;
 202     }
 203   }
 204   if (carry != 0) {
 205     AddCarry(that.digits_, carry);
 206   }
 207   return DivideBy<2>() != 0;
 208 }
 209
 210 template <int PREC, int LOG10RADIX>
 211 void BigRadixFloatingPointNumber<PREC, LOG10RADIX>::Minimize(
 212     BigRadixFloatingPointNumber &&less, BigRadixFloatingPointNumber &&more) {
 213   int leastExponent{exponent_};
 214   if (less.exponent_ < leastExponent) {
 215     leastExponent = less.exponent_;
 216   }
 217   if (more.exponent_ < leastExponent) {
 218     leastExponent = more.exponent_;
 219   }
 220   while (exponent_ > leastExponent) {
 221     --exponent_;
 222     MultiplyBy<10>();
 223   }
 224   while (less.exponent_ > leastExponent) {
 225     --less.exponent_;
 226     less.MultiplyBy<10>();
 227   }
 228   while (more.exponent_ > leastExponent) {
 229     --more.exponent_;
 230     more.MultiplyBy<10>();
 231   }
 232   if (less.Mean(*this)) {
 233     less.AddCarry(); // round up
 234   }
 235   if (!more.Mean(*this)) {
 236     more.Decrement(); // round down
 237   }
 238   while (less.digits_ < more.digits_) {
 239     less.digit_[less.digits_++] = 0;
 240   }
 241   while (more.digits_ < less.digits_) {
 242     more.digit_[more.digits_++] = 0;
 243   }
 244   int digits{more.digits_};
 245   int same{0};
 246   while (same < digits &&
 247       less.digit_[digits - 1 - same] == more.digit_[digits - 1 - same]) {
 248     ++same;
 249   }
 250   if (same == digits) {
 251     return;
 252   }
 253   digits_ = same + 1;
 254   int offset{digits - digits_};
 255   exponent_ += offset * log10Radix;
 256   for (int j{0}; j < digits_; ++j) {
 257     digit_[j] = more.digit_[j + offset];
 258   }
 259   Digit least{less.digit_[offset]};
 260   Digit my{digit_[0]};
 261   while (true) {
 262     Digit q{my / 10u};
 263     Digit r{my - 10 * q};
 264     Digit lq{least / 10u};
 265     Digit lr{least - 10 * lq};
 266     if (r != 0 && lq == q) {
 267       Digit sub{(r - lr) >> 1};
 268       digit_[0] -= sub;
 269       break;
 270     } else {
 271       least = lq;
 272       my = q;
 273       DivideBy<10>();
 274       ++exponent_;
 275     }
 276   }
 277   Normalize();
 278 }
 279
 280 template <int PREC>
 281 ConversionToDecimalResult ConvertToDecimal(char *buffer, std::size_t size,
 282     enum DecimalConversionFlags flags, int digits,
 283     enum FortranRounding rounding, BinaryFloatingPointNumber<PREC> x) {
 284   if (x.IsNaN()) {
 285     return {"NaN", 3, 0, Invalid};
 286   } else if (x.IsInfinite()) {
 287     if (x.IsNegative()) {
 288       return {"-Inf", 4, 0, Exact};
 289     } else if (flags & AlwaysSign) {
 290       return {"+Inf", 4, 0, Exact};
 291     } else {
 292       return {"Inf", 3, 0, Exact};
 293     }
 294   } else {
 295     using Big = BigRadixFloatingPointNumber<PREC>;
 296     Big number{x, rounding};
 297     if ((flags & Minimize) && !x.IsZero()) {
 298       // To emit the fewest decimal digits necessary to represent the value
 299       // in such a way that decimal-to-binary conversion to the same format
 300       // with a fixed assumption about rounding will return the same binary
 301       // value, we also perform binary-to-decimal conversion on the two
 302       // binary values immediately adjacent to this one, use them to identify
 303       // the bounds of the range of decimal values that will map back to the
 304       // original binary value, and find a (not necessary unique) shortest
 305       // decimal sequence in that range.
 306       using Binary = typename Big::Real;
 307       Binary less{x};
 308       less.Previous();
 309       Binary more{x};
 310       if (!x.IsMaximalFiniteMagnitude()) {
 311         more.Next();
 312       }
 313       number.Minimize(Big{less, rounding}, Big{more, rounding});
 314     }
 315     return number.ConvertToDecimal(buffer, size, flags, digits);
 316   }
 317 }
 318
 319 template ConversionToDecimalResult ConvertToDecimal<8>(char *, std::size_t,
 320     enum DecimalConversionFlags, int, enum FortranRounding,
 321     BinaryFloatingPointNumber<8>);
 322 template ConversionToDecimalResult ConvertToDecimal<11>(char *, std::size_t,
 323     enum DecimalConversionFlags, int, enum FortranRounding,
 324     BinaryFloatingPointNumber<11>);
 325 template ConversionToDecimalResult ConvertToDecimal<24>(char *, std::size_t,
 326     enum DecimalConversionFlags, int, enum FortranRounding,
 327     BinaryFloatingPointNumber<24>);
 328 template ConversionToDecimalResult ConvertToDecimal<53>(char *, std::size_t,
 329     enum DecimalConversionFlags, int, enum FortranRounding,
 330     BinaryFloatingPointNumber<53>);
 331 template ConversionToDecimalResult ConvertToDecimal<64>(char *, std::size_t,
 332     enum DecimalConversionFlags, int, enum FortranRounding,
 333     BinaryFloatingPointNumber<64>);
 334 template ConversionToDecimalResult ConvertToDecimal<113>(char *, std::size_t,
 335     enum DecimalConversionFlags, int, enum FortranRounding,
 336     BinaryFloatingPointNumber<113>);
 337
 338 extern "C" {
 339 ConversionToDecimalResult ConvertFloatToDecimal(char *buffer, std::size_t size,
 340     enum DecimalConversionFlags flags, int digits,
 341     enum FortranRounding rounding, float x) {
 342   return Fortran::decimal::ConvertToDecimal(buffer, size, flags, digits,
 343       rounding, Fortran::decimal::BinaryFloatingPointNumber<24>(x));
 344 }
 345
 346 ConversionToDecimalResult ConvertDoubleToDecimal(char *buffer, std::size_t size,
 347     enum DecimalConversionFlags flags, int digits,
 348     enum FortranRounding rounding, double x) {
 349   return Fortran::decimal::ConvertToDecimal(buffer, size, flags, digits,
 350       rounding, Fortran::decimal::BinaryFloatingPointNumber<53>(x));
 351 }
 352
 353 #if LDBL_MANT_DIG == 64
 354 ConversionToDecimalResult ConvertLongDoubleToDecimal(char *buffer,
 355     std::size_t size, enum DecimalConversionFlags flags, int digits,
 356     enum FortranRounding rounding, long double x) {
 357   return Fortran::decimal::ConvertToDecimal(buffer, size, flags, digits,
 358       rounding, Fortran::decimal::BinaryFloatingPointNumber<64>(x));
 359 }
 360 #elif LDBL_MANT_DIG == 113
 361 ConversionToDecimalResult ConvertLongDoubleToDecimal(char *buffer,
 362     std::size_t size, enum DecimalConversionFlags flags, int digits,
 363     enum FortranRounding rounding, long double x) {
 364   return Fortran::decimal::ConvertToDecimal(buffer, size, flags, digits,
 365       rounding, Fortran::decimal::BinaryFloatingPointNumber<113>(x));
 366 }
 367 #endif
 368 }
 369
 370 template <int PREC, int LOG10RADIX>
 371 template <typename STREAM>
 372 STREAM &BigRadixFloatingPointNumber<PREC, LOG10RADIX>::Dump(STREAM &o) const {
 373   if (isNegative_) {
 374     o << '-';
 375   }
 376   o << "10**(" << exponent_ << ") * ...\n";
 377   for (int j{digits_}; --j >= 0;) {
 378     std::string str{std::to_string(digit_[j])};
 379     o << std::string(20 - str.size(), ' ') << str << " [" << j << ']';
 380     if (j + 1 == digitLimit_) {
 381       o << " (limit)";
 382     }
 383     o << '\n';
 384   }
 385   return o;
 386 }
 387 } // namespace Fortran::decimal