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