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 <string>
13 
14 namespace Fortran::decimal {
15 
16 template <int PREC, int LOG10RADIX>
17 BigRadixFloatingPointNumber<PREC, LOG10RADIX>::BigRadixFloatingPointNumber(
18     BinaryFloatingPointNumber<PREC> x, enum FortranRounding rounding)
19     : rounding_{rounding} {
20   bool negative{x.IsNegative()};
21   if (x.IsZero()) {
22     isNegative_ = negative;
23     return;
24   }
25   if (negative) {
26     x.Negate();
27   }
28   int twoPow{x.UnbiasedExponent()};
29   twoPow -= x.bits - 1;
30   if (!x.isImplicitMSB) {
31     ++twoPow;
32   }
33   int lshift{x.exponentBits};
34   if (twoPow <= -lshift) {
35     twoPow += lshift;
36     lshift = 0;
37   } else if (twoPow < 0) {
38     lshift += twoPow;
39     twoPow = 0;
40   }
41   auto word{x.Fraction()};
42   word <<= lshift;
43   SetTo(word);
44   isNegative_ = negative;
45 
46   // The significand is now encoded in *this as an integer (D) and
47   // decimal exponent (E):  x = D * 10.**E * 2.**twoPow
48   // twoPow can be positive or negative.
49   // The goal now is to get twoPow up or down to zero, leaving us with
50   // only decimal digits and decimal exponent.  This is done by
51   // fast multiplications and divisions of D by 2 and 5.
52 
53   // (5*D) * 10.**E * 2.**twoPow -> D * 10.**(E+1) * 2.**(twoPow-1)
54   for (; twoPow > 0 && IsDivisibleBy<5>(); --twoPow) {
55     DivideBy<5>();
56     ++exponent_;
57   }
58 
59   int overflow{0};
60   for (; twoPow >= 9; twoPow -= 9) {
61     // D * 10.**E * 2.**twoPow -> (D*(2**9)) * 10.**E * 2.**(twoPow-9)
62     overflow |= MultiplyBy<512>();
63   }
64   for (; twoPow >= 3; twoPow -= 3) {
65     // D * 10.**E * 2.**twoPow -> (D*(2**3)) * 10.**E * 2.**(twoPow-3)
66     overflow |= MultiplyBy<8>();
67   }
68   for (; twoPow > 0; --twoPow) {
69     // D * 10.**E * 2.**twoPow -> (2*D) * 10.**E * 2.**(twoPow-1)
70     overflow |= MultiplyBy<2>();
71   }
72 
73   overflow |= DivideByPowerOfTwoInPlace(-twoPow);
74   assert(overflow == 0);
75   Normalize();
76 }
77 
78 template <int PREC, int LOG10RADIX>
79 ConversionToDecimalResult
80 BigRadixFloatingPointNumber<PREC, LOG10RADIX>::ConvertToDecimal(char *buffer,
81     std::size_t n, enum DecimalConversionFlags flags, int maxDigits) const {
82   if (n < static_cast<std::size_t>(3 + digits_ * LOG10RADIX)) {
83     return {nullptr, 0, 0, Overflow};
84   }
85   char *start{buffer};
86   if (isNegative_) {
87     *start++ = '-';
88   } else if (flags & AlwaysSign) {
89     *start++ = '+';
90   }
91   if (IsZero()) {
92     *start++ = '0';
93     *start = '\0';
94     return {buffer, static_cast<std::size_t>(start - buffer), 0, Exact};
95   }
96   char *p{start};
97   static_assert((LOG10RADIX % 2) == 0, "radix not a power of 100");
98   static const char lut[] = "0001020304050607080910111213141516171819"
99                             "2021222324252627282930313233343536373839"
100                             "4041424344454647484950515253545556575859"
101                             "6061626364656667686970717273747576777879"
102                             "8081828384858687888990919293949596979899";
103   // Treat the MSD specially: don't emit leading zeroes.
104   Digit dig{digit_[digits_ - 1]};
105   char stack[LOG10RADIX], *sp{stack};
106   for (int k{0}; k < log10Radix; k += 2) {
107     Digit newDig{dig / 100};
108     auto d{static_cast<std::uint32_t>(dig) -
109         std::uint32_t{100} * static_cast<std::uint32_t>(newDig)};
110     dig = newDig;
111     const char *q{lut + d + d};
112     *sp++ = q[1];
113     *sp++ = q[0];
114   }
115   while (sp > stack && sp[-1] == '0') {
116     --sp;
117   }
118   while (sp > stack) {
119     *p++ = *--sp;
120   }
121   for (int j{digits_ - 1}; j-- > 0;) {
122     Digit dig{digit_[j]};
123     char *reverse{p += log10Radix};
124     for (int k{0}; k < log10Radix; k += 2) {
125       Digit newDig{dig / 100};
126       auto d{static_cast<std::uint32_t>(dig) -
127           std::uint32_t{100} * static_cast<std::uint32_t>(newDig)};
128       dig = newDig;
129       const char *q{lut + d + d};
130       *--reverse = q[1];
131       *--reverse = q[0];
132     }
133   }
134   // Adjust exponent so the effective decimal point is to
135   // the left of the first digit.
136   int expo = exponent_ + p - start;
137   // Trim trailing zeroes.
138   while (p[-1] == '0') {
139     --p;
140   }
141   char *end{start + maxDigits};
142   if (maxDigits == 0) {
143     p = end;
144   }
145   if (p <= end) {
146     *p = '\0';
147     return {buffer, static_cast<std::size_t>(p - buffer), expo, Exact};
148   } else {
149     // Apply a digit limit, possibly with rounding.
150     bool incr{false};
151     switch (rounding_) {
152     case RoundNearest:
153       incr = *end > '5' ||
154           (*end == '5' && (p > end + 1 || ((end[-1] - '0') & 1) != 0));
155       break;
156     case RoundUp:
157       incr = !isNegative_;
158       break;
159     case RoundDown:
160       incr = isNegative_;
161       break;
162     case RoundToZero:
163       break;
164     case RoundCompatible:
165       incr = *end >= '5';
166       break;
167     }
168     p = end;
169     if (incr) {
170       while (p > start && p[-1] == '9') {
171         --p;
172       }
173       if (p == start) {
174         *p++ = '1';
175         ++expo;
176       } else {
177         ++p[-1];
178       }
179     }
180 
181     *p = '\0';
182     return {buffer, static_cast<std::size_t>(p - buffer), expo, Inexact};
183   }
184 }
185 
186 template <int PREC, int LOG10RADIX>
187 bool BigRadixFloatingPointNumber<PREC, LOG10RADIX>::Mean(
188     const BigRadixFloatingPointNumber &that) {
189   while (digits_ < that.digits_) {
190     digit_[digits_++] = 0;
191   }
192   int carry{0};
193   for (int j{0}; j < that.digits_; ++j) {
194     Digit v{digit_[j] + that.digit_[j] + carry};
195     if (v >= radix) {
196       digit_[j] = v - radix;
197       carry = 1;
198     } else {
199       digit_[j] = v;
200       carry = 0;
201     }
202   }
203   if (carry != 0) {
204     AddCarry(that.digits_, carry);
205   }
206   return DivideBy<2>() != 0;
207 }
208 
209 template <int PREC, int LOG10RADIX>
210 void BigRadixFloatingPointNumber<PREC, LOG10RADIX>::Minimize(
211     BigRadixFloatingPointNumber &&less, BigRadixFloatingPointNumber &&more) {
212   int leastExponent{exponent_};
213   if (less.exponent_ < leastExponent) {
214     leastExponent = less.exponent_;
215   }
216   if (more.exponent_ < leastExponent) {
217     leastExponent = more.exponent_;
218   }
219   while (exponent_ > leastExponent) {
220     --exponent_;
221     MultiplyBy<10>();
222   }
223   while (less.exponent_ > leastExponent) {
224     --less.exponent_;
225     less.MultiplyBy<10>();
226   }
227   while (more.exponent_ > leastExponent) {
228     --more.exponent_;
229     more.MultiplyBy<10>();
230   }
231   if (less.Mean(*this)) {
232     less.AddCarry(); // round up
233   }
234   if (!more.Mean(*this)) {
235     more.Decrement(); // round down
236   }
237   while (less.digits_ < more.digits_) {
238     less.digit_[less.digits_++] = 0;
239   }
240   while (more.digits_ < less.digits_) {
241     more.digit_[more.digits_++] = 0;
242   }
243   int digits{more.digits_};
244   int same{0};
245   while (same < digits &&
246       less.digit_[digits - 1 - same] == more.digit_[digits - 1 - same]) {
247     ++same;
248   }
249   if (same == digits) {
250     return;
251   }
252   digits_ = same + 1;
253   int offset{digits - digits_};
254   exponent_ += offset * log10Radix;
255   for (int j{0}; j < digits_; ++j) {
256     digit_[j] = more.digit_[j + offset];
257   }
258   Digit least{less.digit_[offset]};
259   Digit my{digit_[0]};
260   while (true) {
261     Digit q{my / 10u};
262     Digit r{my - 10 * q};
263     Digit lq{least / 10u};
264     Digit lr{least - 10 * lq};
265     if (r != 0 && lq == q) {
266       Digit sub{(r - lr) >> 1};
267       digit_[0] -= sub;
268       break;
269     } else {
270       least = lq;
271       my = q;
272       DivideBy<10>();
273       ++exponent_;
274     }
275   }
276   Normalize();
277 }
278 
279 template <int PREC>
280 ConversionToDecimalResult ConvertToDecimal(char *buffer, std::size_t size,
281     enum DecimalConversionFlags flags, int digits,
282     enum FortranRounding rounding, BinaryFloatingPointNumber<PREC> x) {
283   if (x.IsNaN()) {
284     return {"NaN", 3, 0, Invalid};
285   } else if (x.IsInfinite()) {
286     if (x.IsNegative()) {
287       return {"-Inf", 4, 0, Exact};
288     } else if (flags & AlwaysSign) {
289       return {"+Inf", 4, 0, Exact};
290     } else {
291       return {"Inf", 3, 0, Exact};
292     }
293   } else {
294     using Big = BigRadixFloatingPointNumber<PREC>;
295     Big number{x, rounding};
296     if ((flags & Minimize) && !x.IsZero()) {
297       // To emit the fewest decimal digits necessary to represent the value
298       // in such a way that decimal-to-binary conversion to the same format
299       // with a fixed assumption about rounding will return the same binary
300       // value, we also perform binary-to-decimal conversion on the two
301       // binary values immediately adjacent to this one, use them to identify
302       // the bounds of the range of decimal values that will map back to the
303       // original binary value, and find a (not necessary unique) shortest
304       // decimal sequence in that range.
305       using Binary = typename Big::Real;
306       Binary less{x};
307       less.Previous();
308       Binary more{x};
309       if (!x.IsMaximalFiniteMagnitude()) {
310         more.Next();
311       }
312       number.Minimize(Big{less, rounding}, Big{more, rounding});
313     }
314     return number.ConvertToDecimal(buffer, size, flags, digits);
315   }
316 }
317 
318 template ConversionToDecimalResult ConvertToDecimal<8>(char *, std::size_t,
319     enum DecimalConversionFlags, int, enum FortranRounding,
320     BinaryFloatingPointNumber<8>);
321 template ConversionToDecimalResult ConvertToDecimal<11>(char *, std::size_t,
322     enum DecimalConversionFlags, int, enum FortranRounding,
323     BinaryFloatingPointNumber<11>);
324 template ConversionToDecimalResult ConvertToDecimal<24>(char *, std::size_t,
325     enum DecimalConversionFlags, int, enum FortranRounding,
326     BinaryFloatingPointNumber<24>);
327 template ConversionToDecimalResult ConvertToDecimal<53>(char *, std::size_t,
328     enum DecimalConversionFlags, int, enum FortranRounding,
329     BinaryFloatingPointNumber<53>);
330 template ConversionToDecimalResult ConvertToDecimal<64>(char *, std::size_t,
331     enum DecimalConversionFlags, int, enum FortranRounding,
332     BinaryFloatingPointNumber<64>);
333 template ConversionToDecimalResult ConvertToDecimal<113>(char *, std::size_t,
334     enum DecimalConversionFlags, int, enum FortranRounding,
335     BinaryFloatingPointNumber<113>);
336 
337 extern "C" {
338 ConversionToDecimalResult ConvertFloatToDecimal(char *buffer, std::size_t size,
339     enum DecimalConversionFlags flags, int digits,
340     enum FortranRounding rounding, float x) {
341   return Fortran::decimal::ConvertToDecimal(buffer, size, flags, digits,
342       rounding, Fortran::decimal::BinaryFloatingPointNumber<24>(x));
343 }
344 
345 ConversionToDecimalResult ConvertDoubleToDecimal(char *buffer, std::size_t size,
346     enum DecimalConversionFlags flags, int digits,
347     enum FortranRounding rounding, double x) {
348   return Fortran::decimal::ConvertToDecimal(buffer, size, flags, digits,
349       rounding, Fortran::decimal::BinaryFloatingPointNumber<53>(x));
350 }
351 
352 #if LONG_DOUBLE == 80
353 ConversionToDecimalResult ConvertLongDoubleToDecimal(char *buffer,
354     std::size_t size, enum DecimalConversionFlags flags, int digits,
355     enum FortranRounding rounding, long double x) {
356   return Fortran::decimal::ConvertToDecimal(buffer, size, flags, digits,
357       rounding, Fortran::decimal::BinaryFloatingPointNumber<64>(x));
358 }
359 #elif LONG_DOUBLE == 128
360 ConversionToDecimalResult ConvertLongDoubleToDecimal(char *buffer,
361     std::size_t size, enum DecimalConversionFlags flags, int digits,
362     enum FortranRounding rounding, long double x) {
363   return Fortran::decimal::ConvertToDecimal(buffer, size, flags, digits,
364       rounding, Fortran::decimal::BinaryFloatingPointNumber<113>(x));
365 }
366 #endif
367 }
368 
369 template <int PREC, int LOG10RADIX>
370 template <typename STREAM>
371 STREAM &BigRadixFloatingPointNumber<PREC, LOG10RADIX>::Dump(STREAM &o) const {
372   if (isNegative_) {
373     o << '-';
374   }
375   o << "10**(" << exponent_ << ") * ...\n";
376   for (int j{digits_}; --j >= 0;) {
377     std::string str{std::to_string(digit_[j])};
378     o << std::string(20 - str.size(), ' ') << str << " [" << j << ']';
379     if (j + 1 == digitLimit_) {
380       o << " (limit)";
381     }
382     o << '\n';
383   }
384   return o;
385 }
386 } // namespace Fortran::decimal
387