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