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