393 lines · 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-exception6//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 15namespace Fortran::decimal {16 17template <int PREC, int LOG10RADIX>18BigRadixFloatingPointNumber<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) and48 // decimal exponent (E): x = D * 10.**E * 2.**twoPow49 // twoPow can be positive or negative.50 // The goal now is to get twoPow up or down to zero, leaving us with51 // only decimal digits and decimal exponent. This is done by52 // 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 79template <int PREC, int LOG10RADIX>80ConversionToDecimalResult81BigRadixFloatingPointNumber<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 to136 // 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 187template <int PREC, int LOG10RADIX>188bool 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 210template <int PREC, int LOG10RADIX>211void 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 up234 }235 if (!more.Mean(*this)) {236 more.Decrement(); // round down237 }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 280template <int PREC>281ConversionToDecimalResult 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 value299 // in such a way that decimal-to-binary conversion to the same format300 // with a fixed assumption about rounding will return the same binary301 // value, we also perform binary-to-decimal conversion on the two302 // binary values immediately adjacent to this one, use them to identify303 // the bounds of the range of decimal values that will map back to the304 // original binary value, and find a (not necessary unique) shortest305 // 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 319template ConversionToDecimalResult ConvertToDecimal<8>(char *, std::size_t,320 enum DecimalConversionFlags, int, enum FortranRounding,321 BinaryFloatingPointNumber<8>);322template ConversionToDecimalResult ConvertToDecimal<11>(char *, std::size_t,323 enum DecimalConversionFlags, int, enum FortranRounding,324 BinaryFloatingPointNumber<11>);325template ConversionToDecimalResult ConvertToDecimal<24>(char *, std::size_t,326 enum DecimalConversionFlags, int, enum FortranRounding,327 BinaryFloatingPointNumber<24>);328template ConversionToDecimalResult ConvertToDecimal<53>(char *, std::size_t,329 enum DecimalConversionFlags, int, enum FortranRounding,330 BinaryFloatingPointNumber<53>);331template ConversionToDecimalResult ConvertToDecimal<64>(char *, std::size_t,332 enum DecimalConversionFlags, int, enum FortranRounding,333 BinaryFloatingPointNumber<64>);334template ConversionToDecimalResult ConvertToDecimal<113>(char *, std::size_t,335 enum DecimalConversionFlags, int, enum FortranRounding,336 BinaryFloatingPointNumber<113>);337 338extern "C" {339RT_EXT_API_GROUP_BEGIN340 341ConversionToDecimalResult ConvertFloatToDecimal(char *buffer, std::size_t size,342 enum DecimalConversionFlags flags, int digits,343 enum FortranRounding rounding, float x) {344 return Fortran::decimal::ConvertToDecimal(buffer, size, flags, digits,345 rounding, Fortran::decimal::BinaryFloatingPointNumber<24>(x));346}347 348ConversionToDecimalResult ConvertDoubleToDecimal(char *buffer, std::size_t size,349 enum DecimalConversionFlags flags, int digits,350 enum FortranRounding rounding, double x) {351 return Fortran::decimal::ConvertToDecimal(buffer, size, flags, digits,352 rounding, Fortran::decimal::BinaryFloatingPointNumber<53>(x));353}354 355#if LDBL_MANT_DIG == 64356ConversionToDecimalResult ConvertLongDoubleToDecimal(char *buffer,357 std::size_t size, enum DecimalConversionFlags flags, int digits,358 enum FortranRounding rounding, long double x) {359 return Fortran::decimal::ConvertToDecimal(buffer, size, flags, digits,360 rounding, Fortran::decimal::BinaryFloatingPointNumber<64>(x));361}362#elif LDBL_MANT_DIG == 113363ConversionToDecimalResult ConvertLongDoubleToDecimal(char *buffer,364 std::size_t size, enum DecimalConversionFlags flags, int digits,365 enum FortranRounding rounding, long double x) {366 return Fortran::decimal::ConvertToDecimal(buffer, size, flags, digits,367 rounding, Fortran::decimal::BinaryFloatingPointNumber<113>(x));368}369#endif370 371RT_EXT_API_GROUP_END372} // extern "C"373 374template <int PREC, int LOG10RADIX>375template <typename STREAM>376STREAM &BigRadixFloatingPointNumber<PREC, LOG10RADIX>::Dump(STREAM &o) const {377 if (isNegative_) {378 o << '-';379 }380 o << "10**(" << exponent_ << ") * ... (rounding "381 << static_cast<int>(rounding_) << ")\n";382 for (int j{digits_}; --j >= 0;) {383 std::string str{std::to_string(digit_[j])};384 o << std::string(20 - str.size(), ' ') << str << " [" << j << ']';385 if (j + 1 == digitLimit_) {386 o << " (limit)";387 }388 o << '\n';389 }390 return o;391}392} // namespace Fortran::decimal393