• Home
  • Line#
  • Scopes#
  • Navigate#
  • Raw
  • Download
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