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   static constexpr Digit hundredth{radix / 100};
104   // Treat the MSD specially: don't emit leading zeroes.
105   Digit dig{digit_[digits_ - 1]};
106   for (int k{0}; k < LOG10RADIX; k += 2) {
107     Digit d{common::DivideUnsignedBy<Digit, hundredth>(dig)};
108     dig = 100 * (dig - d * hundredth);
109     const char *q{lut + 2 * d};
110     if (q[0] != '0' || p > start) {
111       *p++ = q[0];
112       *p++ = q[1];
113     } else if (q[1] != '0') {
114       *p++ = q[1];
115     }
116   }
117   for (int j{digits_ - 1}; j-- > 0;) {
118     Digit dig{digit_[j]};
119     for (int k{0}; k < log10Radix; k += 2) {
120       Digit d{common::DivideUnsignedBy<Digit, hundredth>(dig)};
121       dig = 100 * (dig - d * hundredth);
122       const char *q{lut + 2 * d};
123       *p++ = q[0];
124       *p++ = q[1];
125     }
126   }
127   // Adjust exponent so the effective decimal point is to
128   // the left of the first digit.
129   int expo = exponent_ + p - start;
130   // Trim trailing zeroes.
131   while (p[-1] == '0') {
132     --p;
133   }
134   char *end{start + maxDigits};
135   if (maxDigits == 0) {
136     p = end;
137   }
138   if (p <= end) {
139     *p = '\0';
140     return {buffer, static_cast<std::size_t>(p - buffer), expo, Exact};
141   } else {
142     // Apply a digit limit, possibly with rounding.
143     bool incr{false};
144     switch (rounding_) {
145     case RoundNearest:
146       incr = *end > '5' ||
147           (*end == '5' && (p > end + 1 || ((end[-1] - '0') & 1) != 0));
148       break;
149     case RoundUp:
150       incr = !isNegative_;
151       break;
152     case RoundDown:
153       incr = isNegative_;
154       break;
155     case RoundToZero:
156       break;
157     case RoundCompatible:
158       incr = *end >= '5';
159       break;
160     }
161     p = end;
162     if (incr) {
163       while (p > start && p[-1] == '9') {
164         --p;
165       }
166       if (p == start) {
167         *p++ = '1';
168         ++expo;
169       } else {
170         ++p[-1];
171       }
172     }
173 
174     *p = '\0';
175     return {buffer, static_cast<std::size_t>(p - buffer), expo, Inexact};
176   }
177 }
178 
179 template <int PREC, int LOG10RADIX>
180 bool BigRadixFloatingPointNumber<PREC, LOG10RADIX>::Mean(
181     const BigRadixFloatingPointNumber &that) {
182   while (digits_ < that.digits_) {
183     digit_[digits_++] = 0;
184   }
185   int carry{0};
186   for (int j{0}; j < that.digits_; ++j) {
187     Digit v{digit_[j] + that.digit_[j] + carry};
188     if (v >= radix) {
189       digit_[j] = v - radix;
190       carry = 1;
191     } else {
192       digit_[j] = v;
193       carry = 0;
194     }
195   }
196   if (carry != 0) {
197     AddCarry(that.digits_, carry);
198   }
199   return DivideBy<2>() != 0;
200 }
201 
202 template <int PREC, int LOG10RADIX>
203 void BigRadixFloatingPointNumber<PREC, LOG10RADIX>::Minimize(
204     BigRadixFloatingPointNumber &&less, BigRadixFloatingPointNumber &&more) {
205   int leastExponent{exponent_};
206   if (less.exponent_ < leastExponent) {
207     leastExponent = less.exponent_;
208   }
209   if (more.exponent_ < leastExponent) {
210     leastExponent = more.exponent_;
211   }
212   while (exponent_ > leastExponent) {
213     --exponent_;
214     MultiplyBy<10>();
215   }
216   while (less.exponent_ > leastExponent) {
217     --less.exponent_;
218     less.MultiplyBy<10>();
219   }
220   while (more.exponent_ > leastExponent) {
221     --more.exponent_;
222     more.MultiplyBy<10>();
223   }
224   if (less.Mean(*this)) {
225     less.AddCarry(); // round up
226   }
227   if (!more.Mean(*this)) {
228     more.Decrement(); // round down
229   }
230   while (less.digits_ < more.digits_) {
231     less.digit_[less.digits_++] = 0;
232   }
233   while (more.digits_ < less.digits_) {
234     more.digit_[more.digits_++] = 0;
235   }
236   int digits{more.digits_};
237   int same{0};
238   while (same < digits &&
239       less.digit_[digits - 1 - same] == more.digit_[digits - 1 - same]) {
240     ++same;
241   }
242   if (same == digits) {
243     return;
244   }
245   digits_ = same + 1;
246   int offset{digits - digits_};
247   exponent_ += offset * log10Radix;
248   for (int j{0}; j < digits_; ++j) {
249     digit_[j] = more.digit_[j + offset];
250   }
251   Digit least{less.digit_[offset]};
252   Digit my{digit_[0]};
253   while (true) {
254     Digit q{common::DivideUnsignedBy<Digit, 10>(my)};
255     Digit r{my - 10 * q};
256     Digit lq{common::DivideUnsignedBy<Digit, 10>(least)};
257     Digit lr{least - 10 * lq};
258     if (r != 0 && lq == q) {
259       Digit sub{(r - lr) >> 1};
260       digit_[0] -= sub;
261       break;
262     } else {
263       least = lq;
264       my = q;
265       DivideBy<10>();
266       ++exponent_;
267     }
268   }
269   Normalize();
270 }
271 
272 template <int PREC>
273 ConversionToDecimalResult ConvertToDecimal(char *buffer, std::size_t size,
274     enum DecimalConversionFlags flags, int digits,
275     enum FortranRounding rounding, BinaryFloatingPointNumber<PREC> x) {
276   if (x.IsNaN()) {
277     return {"NaN", 3, 0, Invalid};
278   } else if (x.IsInfinite()) {
279     if (x.IsNegative()) {
280       return {"-Inf", 4, 0, Exact};
281     } else if (flags & AlwaysSign) {
282       return {"+Inf", 4, 0, Exact};
283     } else {
284       return {"Inf", 3, 0, Exact};
285     }
286   } else {
287     using Big = BigRadixFloatingPointNumber<PREC>;
288     Big number{x, rounding};
289     if ((flags & Minimize) && !x.IsZero()) {
290       // To emit the fewest decimal digits necessary to represent the value
291       // in such a way that decimal-to-binary conversion to the same format
292       // with a fixed assumption about rounding will return the same binary
293       // value, we also perform binary-to-decimal conversion on the two
294       // binary values immediately adjacent to this one, use them to identify
295       // the bounds of the range of decimal values that will map back to the
296       // original binary value, and find a (not necessary unique) shortest
297       // decimal sequence in that range.
298       using Binary = typename Big::Real;
299       Binary less{x};
300       less.Previous();
301       Binary more{x};
302       if (!x.IsMaximalFiniteMagnitude()) {
303         more.Next();
304       }
305       number.Minimize(Big{less, rounding}, Big{more, rounding});
306     } else {
307     }
308     return number.ConvertToDecimal(buffer, size, flags, digits);
309   }
310 }
311 
312 template ConversionToDecimalResult ConvertToDecimal<8>(char *, std::size_t,
313     enum DecimalConversionFlags, int, enum FortranRounding,
314     BinaryFloatingPointNumber<8>);
315 template ConversionToDecimalResult ConvertToDecimal<11>(char *, std::size_t,
316     enum DecimalConversionFlags, int, enum FortranRounding,
317     BinaryFloatingPointNumber<11>);
318 template ConversionToDecimalResult ConvertToDecimal<24>(char *, std::size_t,
319     enum DecimalConversionFlags, int, enum FortranRounding,
320     BinaryFloatingPointNumber<24>);
321 template ConversionToDecimalResult ConvertToDecimal<53>(char *, std::size_t,
322     enum DecimalConversionFlags, int, enum FortranRounding,
323     BinaryFloatingPointNumber<53>);
324 template ConversionToDecimalResult ConvertToDecimal<64>(char *, std::size_t,
325     enum DecimalConversionFlags, int, enum FortranRounding,
326     BinaryFloatingPointNumber<64>);
327 template ConversionToDecimalResult ConvertToDecimal<113>(char *, std::size_t,
328     enum DecimalConversionFlags, int, enum FortranRounding,
329     BinaryFloatingPointNumber<113>);
330 
331 extern "C" {
332 ConversionToDecimalResult ConvertFloatToDecimal(char *buffer, std::size_t size,
333     enum DecimalConversionFlags flags, int digits,
334     enum FortranRounding rounding, float x) {
335   return Fortran::decimal::ConvertToDecimal(buffer, size, flags, digits,
336       rounding, Fortran::decimal::BinaryFloatingPointNumber<24>(x));
337 }
338 
339 ConversionToDecimalResult ConvertDoubleToDecimal(char *buffer, std::size_t size,
340     enum DecimalConversionFlags flags, int digits,
341     enum FortranRounding rounding, double x) {
342   return Fortran::decimal::ConvertToDecimal(buffer, size, flags, digits,
343       rounding, Fortran::decimal::BinaryFloatingPointNumber<53>(x));
344 }
345 
346 #if __x86_64__ && !defined(_MSC_VER)
347 ConversionToDecimalResult ConvertLongDoubleToDecimal(char *buffer,
348     std::size_t size, enum DecimalConversionFlags flags, int digits,
349     enum FortranRounding rounding, long double x) {
350   return Fortran::decimal::ConvertToDecimal(buffer, size, flags, digits,
351       rounding, Fortran::decimal::BinaryFloatingPointNumber<64>(x));
352 }
353 #endif
354 }
355 
356 template <int PREC, int LOG10RADIX>
357 template <typename STREAM>
358 STREAM &BigRadixFloatingPointNumber<PREC, LOG10RADIX>::Dump(STREAM &o) const {
359   if (isNegative_) {
360     o << '-';
361   }
362   o << "10**(" << exponent_ << ") * ...\n";
363   for (int j{digits_}; --j >= 0;) {
364     std::string str{std::to_string(digit_[j])};
365     o << std::string(20 - str.size(), ' ') << str << " [" << j << ']';
366     if (j + 1 == digitLimit_) {
367       o << " (limit)";
368     }
369     o << '\n';
370   }
371   return o;
372 }
373 } // namespace Fortran::decimal
374