| 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 | RT_EXT_API_GROUP_BEGIN |
| 340 | |
| 341 | ConversionToDecimalResult 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 | |
| 348 | ConversionToDecimalResult 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 == 64 |
| 356 | ConversionToDecimalResult 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 == 113 |
| 363 | ConversionToDecimalResult 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 | #endif |
| 370 | |
| 371 | RT_EXT_API_GROUP_END |
| 372 | } // extern "C" |
| 373 | |
| 374 | template <int PREC, int LOG10RADIX> |
| 375 | template <typename STREAM> |
| 376 | STREAM &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::decimal |
| 393 | |