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