397 lines · c
1//===-- lib/Decimal/big-radix-floating-point.h ------------------*- C++ -*-===//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#ifndef FORTRAN_DECIMAL_BIG_RADIX_FLOATING_POINT_H_10#define FORTRAN_DECIMAL_BIG_RADIX_FLOATING_POINT_H_11 12// This is a helper class for use in floating-point conversions between13// binary and decimal representations. It holds a multiple-precision14// integer value using digits of a radix that is a large even power of ten15// (10,000,000,000,000,000 by default, 10**16). These digits are accompanied16// by a signed exponent that denotes multiplication by a power of ten.17// The effective radix point is to the right of the digits (i.e., they do18// not represent a fraction).19//20// The operations supported by this class are limited to those required21// for conversions between binary and decimal representations; it is not22// a general-purpose facility.23 24#include "flang/Common/bit-population-count.h"25#include "flang/Common/leading-zero-bit-count.h"26#include "flang/Common/uint128.h"27#include "flang/Decimal/binary-floating-point.h"28#include "flang/Decimal/decimal.h"29#include <cinttypes>30#include <limits>31#include <type_traits>32 33// Some environments, viz. glibc 2.17, allow the macro HUGE34// to leak out of <math.h>.35#undef HUGE36 37namespace Fortran::decimal {38 39static constexpr std::uint64_t TenToThe(int power) {40 return power <= 0 ? 1 : 10 * TenToThe(power - 1);41}42 43// 10**(LOG10RADIX + 3) must be < 2**wordbits, and LOG10RADIX must be44// even, so that pairs of decimal digits do not straddle Digits.45// So LOG10RADIX must be 16 or 6.46template <int PREC, int LOG10RADIX = 16> class BigRadixFloatingPointNumber {47public:48 using Real = BinaryFloatingPointNumber<PREC>;49 static constexpr int log10Radix{LOG10RADIX};50 51private:52 static constexpr std::uint64_t uint64Radix{TenToThe(log10Radix)};53 static constexpr int minDigitBits{54 64 - common::LeadingZeroBitCount(uint64Radix)};55 using Digit = common::HostUnsignedIntType<minDigitBits>;56 static constexpr Digit radix{uint64Radix};57 static_assert(radix < std::numeric_limits<Digit>::max() / 1000,58 "radix is somehow too big");59 static_assert(radix > std::numeric_limits<Digit>::max() / 10000,60 "radix is somehow too small");61 62 // The base-2 logarithm of the least significant bit that can arise63 // in a subnormal IEEE floating-point number.64 static constexpr int minLog2AnyBit{65 -Real::exponentBias - Real::binaryPrecision};66 67 // The number of Digits needed to represent the smallest subnormal.68 static constexpr int maxDigits{3 - minLog2AnyBit / log10Radix};69 70public:71 explicit RT_API_ATTRS BigRadixFloatingPointNumber(72 enum FortranRounding rounding = RoundNearest)73 : rounding_{rounding} {}74 75 // Converts a binary floating point value.76 explicit RT_API_ATTRS BigRadixFloatingPointNumber(77 Real, enum FortranRounding = RoundNearest);78 79 RT_API_ATTRS BigRadixFloatingPointNumber &SetToZero() {80 isNegative_ = false;81 digits_ = 0;82 exponent_ = 0;83 return *this;84 }85 86 RT_API_ATTRS bool IsInteger() const { return exponent_ >= 0; }87 88 // Converts decimal floating-point to binary.89 RT_API_ATTRS ConversionToBinaryResult<PREC> ConvertToBinary();90 91 // Parses and converts to binary. Handles leading spaces,92 // "NaN", & optionally-signed "Inf". Does not skip internal93 // spaces.94 // The argument is a reference to a pointer that is left95 // pointing to the first character that wasn't parsed.96 RT_API_ATTRS ConversionToBinaryResult<PREC> ConvertToBinary(97 const char *&, const char *end = nullptr);98 99 // Formats a decimal floating-point number to a user buffer.100 // May emit "NaN" or "Inf", or an possibly-signed integer.101 // No decimal point is written, but if it were, it would be102 // after the last digit; the effective decimal exponent is103 // returned as part of the result structure so that it can be104 // formatted by the client.105 RT_API_ATTRS ConversionToDecimalResult ConvertToDecimal(106 char *, std::size_t, enum DecimalConversionFlags, int digits) const;107 108 // Discard decimal digits not needed to distinguish this value109 // from the decimal encodings of two others (viz., the nearest binary110 // floating-point numbers immediately below and above this one).111 // The last decimal digit may not be uniquely determined in all112 // cases, and will be the mean value when that is so (e.g., if113 // last decimal digit values 6-8 would all work, it'll be a 7).114 // This minimization necessarily assumes that the value will be115 // emitted and read back into the same (or less precise) format116 // with default rounding to the nearest value.117 RT_API_ATTRS void Minimize(118 BigRadixFloatingPointNumber &&less, BigRadixFloatingPointNumber &&more);119 120 template <typename STREAM> STREAM &Dump(STREAM &) const;121 122private:123 RT_API_ATTRS BigRadixFloatingPointNumber(124 const BigRadixFloatingPointNumber &that)125 : digits_{that.digits_}, exponent_{that.exponent_},126 isNegative_{that.isNegative_}, rounding_{that.rounding_} {127 for (int j{0}; j < digits_; ++j) {128 digit_[j] = that.digit_[j];129 }130 }131 132 RT_API_ATTRS bool IsZero() const {133 // Don't assume normalization.134 for (int j{0}; j < digits_; ++j) {135 if (digit_[j] != 0) {136 return false;137 }138 }139 return true;140 }141 142 // Predicate: true when 10*value would cause a carry.143 // (When this happens during decimal-to-binary conversion,144 // there are more digits in the input string than can be145 // represented precisely.)146 RT_API_ATTRS bool IsFull() const {147 return digits_ == digitLimit_ && digit_[digits_ - 1] >= radix / 10;148 }149 150 // Sets *this to an unsigned integer value.151 // Returns any remainder.152 template <typename UINT> RT_API_ATTRS UINT SetTo(UINT n) {153 static_assert(154 std::is_same_v<UINT, common::uint128_t> || std::is_unsigned_v<UINT>);155 SetToZero();156 while (n != 0) {157 auto q{n / 10u};158 if (n != q * 10) {159 break;160 }161 ++exponent_;162 n = q;163 }164 if constexpr (sizeof n < sizeof(Digit)) {165 if (n != 0) {166 digit_[digits_++] = n;167 }168 return 0;169 } else {170 while (n != 0 && digits_ < digitLimit_) {171 auto q{n / radix};172 digit_[digits_++] = static_cast<Digit>(n - q * radix);173 n = q;174 }175 return n;176 }177 }178 179 RT_API_ATTRS int RemoveLeastOrderZeroDigits() {180 int remove{0};181 if (digits_ > 0 && digit_[0] == 0) {182 while (remove < digits_ && digit_[remove] == 0) {183 ++remove;184 }185 if (remove >= digits_) {186 digits_ = 0;187 } else if (remove > 0) {188#if defined __GNUC__ && __GNUC__ < 8189 // (&& j + remove < maxDigits) was added to avoid GCC < 8 build failure190 // on -Werror=array-bounds. This can be removed if -Werror is disable.191 for (int j{0}; j + remove < digits_ && (j + remove < maxDigits); ++j) {192#else193 for (int j{0}; j + remove < digits_; ++j) {194#endif195 digit_[j] = digit_[j + remove];196 }197 digits_ -= remove;198 }199 }200 return remove;201 }202 203 RT_API_ATTRS void RemoveLeadingZeroDigits() {204 while (digits_ > 0 && digit_[digits_ - 1] == 0) {205 --digits_;206 }207 }208 209 RT_API_ATTRS void Normalize() {210 RemoveLeadingZeroDigits();211 exponent_ += RemoveLeastOrderZeroDigits() * log10Radix;212 }213 214 // This limited divisibility test only works for even divisors of the radix,215 // which is fine since it's only ever used with 2 and 5.216 template <int N> RT_API_ATTRS bool IsDivisibleBy() const {217 static_assert(N > 1 && radix % N == 0, "bad modulus");218 return digits_ == 0 || (digit_[0] % N) == 0;219 }220 221 template <unsigned DIVISOR> RT_API_ATTRS int DivideBy() {222 Digit remainder{0};223 for (int j{digits_ - 1}; j >= 0; --j) {224 Digit q{digit_[j] / DIVISOR};225 Digit nrem{digit_[j] - DIVISOR * q};226 digit_[j] = q + (radix / DIVISOR) * remainder;227 remainder = nrem;228 }229 return remainder;230 }231 232 RT_API_ATTRS void DivideByPowerOfTwo(int twoPow) { // twoPow <= log10Radix233 Digit remainder{0};234 auto mask{(Digit{1} << twoPow) - 1};235 auto coeff{radix >> twoPow};236 for (int j{digits_ - 1}; j >= 0; --j) {237 auto nrem{digit_[j] & mask};238 digit_[j] = (digit_[j] >> twoPow) + coeff * remainder;239 remainder = nrem;240 }241 }242 243 // Returns true on overflow244 RT_API_ATTRS bool DivideByPowerOfTwoInPlace(int twoPow) {245 if (digits_ > 0) {246 while (twoPow > 0) {247 int chunk{twoPow > log10Radix ? log10Radix : twoPow};248 if ((digit_[0] & ((Digit{1} << chunk) - 1)) == 0) {249 DivideByPowerOfTwo(chunk);250 twoPow -= chunk;251 continue;252 }253 twoPow -= chunk;254 if (digit_[digits_ - 1] >> chunk != 0) {255 if (digits_ == digitLimit_) {256 return true; // overflow257 }258 digit_[digits_++] = 0;259 }260 auto remainder{digit_[digits_ - 1]};261 exponent_ -= log10Radix;262 auto coeff{radix >> chunk}; // precise; radix is (5*2)**log10Radix263 auto mask{(Digit{1} << chunk) - 1};264 for (int j{digits_ - 1}; j >= 1; --j) {265 digit_[j] = (digit_[j - 1] >> chunk) + coeff * remainder;266 remainder = digit_[j - 1] & mask;267 }268 digit_[0] = coeff * remainder;269 }270 }271 return false; // no overflow272 }273 274 RT_API_ATTRS int AddCarry(int position = 0, int carry = 1) {275 for (; position < digits_; ++position) {276 Digit v{digit_[position] + carry};277 if (v < radix) {278 digit_[position] = v;279 return 0;280 }281 digit_[position] = v - radix;282 carry = 1;283 }284 if (digits_ < digitLimit_) {285 digit_[digits_++] = carry;286 return 0;287 }288 Normalize();289 if (digits_ < digitLimit_) {290 digit_[digits_++] = carry;291 return 0;292 }293 return carry;294 }295 296 RT_API_ATTRS void Decrement() {297 for (int j{0}; digit_[j]-- == 0; ++j) {298 digit_[j] = radix - 1;299 }300 }301 302 template <int N> RT_API_ATTRS int MultiplyByHelper(int carry = 0) {303 for (int j{0}; j < digits_; ++j) {304 auto v{N * digit_[j] + carry};305 carry = v / radix;306 digit_[j] = v - carry * radix; // i.e., v % radix307 }308 return carry;309 }310 311 template <int N> RT_API_ATTRS int MultiplyBy(int carry = 0) {312 if (int newCarry{MultiplyByHelper<N>(carry)}) {313 return AddCarry(digits_, newCarry);314 } else {315 return 0;316 }317 }318 319 template <int N> RT_API_ATTRS int MultiplyWithoutNormalization() {320 if (int carry{MultiplyByHelper<N>(0)}) {321 if (digits_ < digitLimit_) {322 digit_[digits_++] = carry;323 return 0;324 } else {325 return carry;326 }327 } else {328 return 0;329 }330 }331 332 RT_API_ATTRS void LoseLeastSignificantDigit(); // with rounding333 334 RT_API_ATTRS void PushCarry(int carry) {335 if (digits_ == maxDigits && RemoveLeastOrderZeroDigits() == 0) {336 LoseLeastSignificantDigit();337 digit_[digits_ - 1] += carry;338 } else {339 digit_[digits_++] = carry;340 }341 }342 343 // Adds another number and then divides by two.344 // Assumes same exponent and sign.345 // Returns true when the result has effectively been rounded down.346 RT_API_ATTRS bool Mean(const BigRadixFloatingPointNumber &);347 348 // Parses a floating-point number; leaves the pointer reference349 // argument pointing at the next character after what was recognized.350 // The "end" argument can be left null if the caller is sure that the351 // string is properly terminated with an addressable character that352 // can't be in a valid floating-point character.353 RT_API_ATTRS bool ParseNumber(const char *&, bool &inexact, const char *end);354 355 using Raw = typename Real::RawType;356 constexpr RT_API_ATTRS Raw SignBit() const {357 return Raw{isNegative_} << (Real::bits - 1);358 }359 constexpr RT_API_ATTRS Raw Infinity() const {360 Raw result{static_cast<Raw>(Real::maxExponent)};361 result <<= Real::significandBits;362 result |= SignBit();363 if constexpr (Real::bits == 80) { // x87364 result |= Raw{1} << 63;365 }366 return result;367 }368 constexpr RT_API_ATTRS Raw NaN(bool isQuiet = true) {369 Raw result{Real::maxExponent};370 result <<= Real::significandBits;371 result |= SignBit();372 if constexpr (Real::bits == 80) { // x87373 result |= Raw{isQuiet ? 3u : 2u} << 62;374 } else {375 Raw quiet{isQuiet ? Raw{2} : Raw{1}};376 quiet <<= Real::significandBits - 2;377 result |= quiet;378 }379 return result;380 }381 constexpr RT_API_ATTRS Raw HUGE() const {382 Raw result{static_cast<Raw>(Real::maxExponent)};383 result <<= Real::significandBits;384 result |= SignBit();385 return result - 1; // decrement exponent, set all significand bits386 }387 388 Digit digit_[maxDigits]; // in little-endian order: digit_[0] is LSD389 int digits_{0}; // # of elements in digit_[] array; zero when zero390 int digitLimit_{maxDigits}; // precision clamp391 int exponent_{0}; // signed power of ten392 bool isNegative_{false};393 enum FortranRounding rounding_ { RoundNearest };394};395} // namespace Fortran::decimal396#endif397