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-exception
6 //
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
13 // between binary decimal representations. It holds a multiple-precision
14 // integer value using digits of a radix that is a large even power of ten
15 // (10,000,000,000,000,000 by default, 10**16). These digits are accompanied
16 // 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 do
18 // not represent a fraction).
19 //
20 // The operations supported by this class are limited to those required
21 // for conversions between binary and decimal representations; it is not
22 // 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 namespace Fortran::decimal {
34
TenToThe(int power)35 static constexpr std::uint64_t TenToThe(int power) {
36 return power <= 0 ? 1 : 10 * TenToThe(power - 1);
37 }
38
39 // 10**(LOG10RADIX + 3) must be < 2**wordbits, and LOG10RADIX must be
40 // even, so that pairs of decimal digits do not straddle Digits.
41 // So LOG10RADIX must be 16 or 6.
42 template <int PREC, int LOG10RADIX = 16> class BigRadixFloatingPointNumber {
43 public:
44 using Real = BinaryFloatingPointNumber<PREC>;
45 static constexpr int log10Radix{LOG10RADIX};
46
47 private:
48 static constexpr std::uint64_t uint64Radix{TenToThe(log10Radix)};
49 static constexpr int minDigitBits{
50 64 - common::LeadingZeroBitCount(uint64Radix)};
51 using Digit = common::HostUnsignedIntType<minDigitBits>;
52 static constexpr Digit radix{uint64Radix};
53 static_assert(radix < std::numeric_limits<Digit>::max() / 1000,
54 "radix is somehow too big");
55 static_assert(radix > std::numeric_limits<Digit>::max() / 10000,
56 "radix is somehow too small");
57
58 // The base-2 logarithm of the least significant bit that can arise
59 // in a subnormal IEEE floating-point number.
60 static constexpr int minLog2AnyBit{
61 -Real::exponentBias - Real::binaryPrecision};
62
63 // The number of Digits needed to represent the smallest subnormal.
64 static constexpr int maxDigits{3 - minLog2AnyBit / log10Radix};
65
66 public:
67 explicit BigRadixFloatingPointNumber(
68 enum FortranRounding rounding = RoundNearest)
69 : rounding_{rounding} {}
70
71 // Converts a binary floating point value.
72 explicit BigRadixFloatingPointNumber(
73 Real, enum FortranRounding = RoundNearest);
74
SetToZero()75 BigRadixFloatingPointNumber &SetToZero() {
76 isNegative_ = false;
77 digits_ = 0;
78 exponent_ = 0;
79 return *this;
80 }
81
82 // Converts decimal floating-point to binary.
83 ConversionToBinaryResult<PREC> ConvertToBinary();
84
85 // Parses and converts to binary. Handles leading spaces,
86 // "NaN", & optionally-signed "Inf". Does not skip internal
87 // spaces.
88 // The argument is a reference to a pointer that is left
89 // pointing to the first character that wasn't parsed.
90 ConversionToBinaryResult<PREC> ConvertToBinary(const char *&);
91
92 // Formats a decimal floating-point number to a user buffer.
93 // May emit "NaN" or "Inf", or an possibly-signed integer.
94 // No decimal point is written, but if it were, it would be
95 // after the last digit; the effective decimal exponent is
96 // returned as part of the result structure so that it can be
97 // formatted by the client.
98 ConversionToDecimalResult ConvertToDecimal(
99 char *, std::size_t, enum DecimalConversionFlags, int digits) const;
100
101 // Discard decimal digits not needed to distinguish this value
102 // from the decimal encodings of two others (viz., the nearest binary
103 // floating-point numbers immediately below and above this one).
104 // The last decimal digit may not be uniquely determined in all
105 // cases, and will be the mean value when that is so (e.g., if
106 // last decimal digit values 6-8 would all work, it'll be a 7).
107 // This minimization necessarily assumes that the value will be
108 // emitted and read back into the same (or less precise) format
109 // with default rounding to the nearest value.
110 void Minimize(
111 BigRadixFloatingPointNumber &&less, BigRadixFloatingPointNumber &&more);
112
113 template <typename STREAM> STREAM &Dump(STREAM &) const;
114
115 private:
BigRadixFloatingPointNumber(const BigRadixFloatingPointNumber & that)116 BigRadixFloatingPointNumber(const BigRadixFloatingPointNumber &that)
117 : digits_{that.digits_}, exponent_{that.exponent_},
118 isNegative_{that.isNegative_}, rounding_{that.rounding_} {
119 for (int j{0}; j < digits_; ++j) {
120 digit_[j] = that.digit_[j];
121 }
122 }
123
IsZero()124 bool IsZero() const {
125 // Don't assume normalization.
126 for (int j{0}; j < digits_; ++j) {
127 if (digit_[j] != 0) {
128 return false;
129 }
130 }
131 return true;
132 }
133
134 // Predicate: true when 10*value would cause a carry.
135 // (When this happens during decimal-to-binary conversion,
136 // there are more digits in the input string than can be
137 // represented precisely.)
IsFull()138 bool IsFull() const {
139 return digits_ == digitLimit_ && digit_[digits_ - 1] >= radix / 10;
140 }
141
142 // Sets *this to an unsigned integer value.
143 // Returns any remainder.
SetTo(UINT n)144 template <typename UINT> UINT SetTo(UINT n) {
145 static_assert(
146 std::is_same_v<UINT, common::uint128_t> || std::is_unsigned_v<UINT>);
147 SetToZero();
148 while (n != 0) {
149 auto q{n / 10u};
150 if (n != q * 10) {
151 break;
152 }
153 ++exponent_;
154 n = q;
155 }
156 if constexpr (sizeof n < sizeof(Digit)) {
157 if (n != 0) {
158 digit_[digits_++] = n;
159 }
160 return 0;
161 } else {
162 while (n != 0 && digits_ < digitLimit_) {
163 auto q{n / radix};
164 digit_[digits_++] = static_cast<Digit>(n - q * radix);
165 n = q;
166 }
167 return n;
168 }
169 }
170
RemoveLeastOrderZeroDigits()171 int RemoveLeastOrderZeroDigits() {
172 int remove{0};
173 if (digits_ > 0 && digit_[0] == 0) {
174 while (remove < digits_ && digit_[remove] == 0) {
175 ++remove;
176 }
177 if (remove >= digits_) {
178 digits_ = 0;
179 } else if (remove > 0) {
180 #if defined __GNUC__ && __GNUC__ < 8
181 // (&& j + remove < maxDigits) was added to avoid GCC < 8 build failure
182 // on -Werror=array-bounds. This can be removed if -Werror is disable.
183 for (int j{0}; j + remove < digits_ && (j + remove < maxDigits); ++j) {
184 #else
185 for (int j{0}; j + remove < digits_; ++j) {
186 #endif
187 digit_[j] = digit_[j + remove];
188 }
189 digits_ -= remove;
190 }
191 }
192 return remove;
193 }
194
195 void RemoveLeadingZeroDigits() {
196 while (digits_ > 0 && digit_[digits_ - 1] == 0) {
197 --digits_;
198 }
199 }
200
201 void Normalize() {
202 RemoveLeadingZeroDigits();
203 exponent_ += RemoveLeastOrderZeroDigits() * log10Radix;
204 }
205
206 // This limited divisibility test only works for even divisors of the radix,
207 // which is fine since it's only ever used with 2 and 5.
208 template <int N> bool IsDivisibleBy() const {
209 static_assert(N > 1 && radix % N == 0, "bad modulus");
210 return digits_ == 0 || (digit_[0] % N) == 0;
211 }
212
213 template <unsigned DIVISOR> int DivideBy() {
214 Digit remainder{0};
215 for (int j{digits_ - 1}; j >= 0; --j) {
216 Digit q{digit_[j] / DIVISOR};
217 Digit nrem{digit_[j] - DIVISOR * q};
218 digit_[j] = q + (radix / DIVISOR) * remainder;
219 remainder = nrem;
220 }
221 return remainder;
222 }
223
224 void DivideByPowerOfTwo(int twoPow) { // twoPow <= log10Radix
225 Digit remainder{0};
226 auto mask{(Digit{1} << twoPow) - 1};
227 auto coeff{radix >> twoPow};
228 for (int j{digits_ - 1}; j >= 0; --j) {
229 auto nrem{digit_[j] & mask};
230 digit_[j] = (digit_[j] >> twoPow) + coeff * remainder;
231 remainder = nrem;
232 }
233 }
234
235 // Returns true on overflow
236 bool DivideByPowerOfTwoInPlace(int twoPow) {
237 if (digits_ > 0) {
238 while (twoPow > 0) {
239 int chunk{twoPow > log10Radix ? log10Radix : twoPow};
240 if ((digit_[0] & ((Digit{1} << chunk) - 1)) == 0) {
241 DivideByPowerOfTwo(chunk);
242 twoPow -= chunk;
243 continue;
244 }
245 twoPow -= chunk;
246 if (digit_[digits_ - 1] >> chunk != 0) {
247 if (digits_ == digitLimit_) {
248 return true; // overflow
249 }
250 digit_[digits_++] = 0;
251 }
252 auto remainder{digit_[digits_ - 1]};
253 exponent_ -= log10Radix;
254 auto coeff{radix >> chunk}; // precise; radix is (5*2)**log10Radix
255 auto mask{(Digit{1} << chunk) - 1};
256 for (int j{digits_ - 1}; j >= 1; --j) {
257 digit_[j] = (digit_[j - 1] >> chunk) + coeff * remainder;
258 remainder = digit_[j - 1] & mask;
259 }
260 digit_[0] = coeff * remainder;
261 }
262 }
263 return false; // no overflow
264 }
265
266 int AddCarry(int position = 0, int carry = 1) {
267 for (; position < digits_; ++position) {
268 Digit v{digit_[position] + carry};
269 if (v < radix) {
270 digit_[position] = v;
271 return 0;
272 }
273 digit_[position] = v - radix;
274 carry = 1;
275 }
276 if (digits_ < digitLimit_) {
277 digit_[digits_++] = carry;
278 return 0;
279 }
280 Normalize();
281 if (digits_ < digitLimit_) {
282 digit_[digits_++] = carry;
283 return 0;
284 }
285 return carry;
286 }
287
288 void Decrement() {
289 for (int j{0}; digit_[j]-- == 0; ++j) {
290 digit_[j] = radix - 1;
291 }
292 }
293
294 template <int N> int MultiplyByHelper(int carry = 0) {
295 for (int j{0}; j < digits_; ++j) {
296 auto v{N * digit_[j] + carry};
297 carry = v / radix;
298 digit_[j] = v - carry * radix; // i.e., v % radix
299 }
300 return carry;
301 }
302
303 template <int N> int MultiplyBy(int carry = 0) {
304 if (int newCarry{MultiplyByHelper<N>(carry)}) {
305 return AddCarry(digits_, newCarry);
306 } else {
307 return 0;
308 }
309 }
310
311 template <int N> int MultiplyWithoutNormalization() {
312 if (int carry{MultiplyByHelper<N>(0)}) {
313 if (digits_ < digitLimit_) {
314 digit_[digits_++] = carry;
315 return 0;
316 } else {
317 return carry;
318 }
319 } else {
320 return 0;
321 }
322 }
323
324 void LoseLeastSignificantDigit(); // with rounding
325
326 void PushCarry(int carry) {
327 if (digits_ == maxDigits && RemoveLeastOrderZeroDigits() == 0) {
328 LoseLeastSignificantDigit();
329 digit_[digits_ - 1] += carry;
330 } else {
331 digit_[digits_++] = carry;
332 }
333 }
334
335 // Adds another number and then divides by two.
336 // Assumes same exponent and sign.
337 // Returns true when the the result has effectively been rounded down.
338 bool Mean(const BigRadixFloatingPointNumber &);
339
340 bool ParseNumber(const char *&, bool &inexact);
341
342 using Raw = typename Real::RawType;
343 constexpr Raw SignBit() const { return Raw{isNegative_} << (Real::bits - 1); }
344 constexpr Raw Infinity() const {
345 return (Raw{Real::maxExponent} << Real::significandBits) | SignBit();
346 }
347 static constexpr Raw NaN() {
348 return (Raw{Real::maxExponent} << Real::significandBits) |
349 (Raw{1} << (Real::significandBits - 2));
350 }
351
352 Digit digit_[maxDigits]; // in little-endian order: digit_[0] is LSD
353 int digits_{0}; // # of elements in digit_[] array; zero when zero
354 int digitLimit_{maxDigits}; // precision clamp
355 int exponent_{0}; // signed power of ten
356 bool isNegative_{false};
357 enum FortranRounding rounding_ { RoundNearest };
358 };
359 } // namespace Fortran::decimal
360 #endif
361