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>
BigRadixFloatingPointNumber(BinaryFloatingPointNumber<PREC> x,enum FortranRounding rounding)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
ConvertToDecimal(char * buffer,std::size_t n,enum DecimalConversionFlags flags,int maxDigits) const81 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>
Mean(const BigRadixFloatingPointNumber & that)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>
Minimize(BigRadixFloatingPointNumber && less,BigRadixFloatingPointNumber && more)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>
ConvertToDecimal(char * buffer,std::size_t size,enum DecimalConversionFlags flags,int digits,enum FortranRounding rounding,BinaryFloatingPointNumber<PREC> x)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 RT_EXT_API_GROUP_BEGIN
340
ConvertFloatToDecimal(char * buffer,std::size_t size,enum DecimalConversionFlags flags,int digits,enum FortranRounding rounding,float x)341 ConversionToDecimalResult ConvertFloatToDecimal(char *buffer, std::size_t size,
342 enum DecimalConversionFlags flags, int digits,
343 enum FortranRounding rounding, float x) {
344 return Fortran::decimal::ConvertToDecimal(buffer, size, flags, digits,
345 rounding, Fortran::decimal::BinaryFloatingPointNumber<24>(x));
346 }
347
ConvertDoubleToDecimal(char * buffer,std::size_t size,enum DecimalConversionFlags flags,int digits,enum FortranRounding rounding,double x)348 ConversionToDecimalResult ConvertDoubleToDecimal(char *buffer, std::size_t size,
349 enum DecimalConversionFlags flags, int digits,
350 enum FortranRounding rounding, double x) {
351 return Fortran::decimal::ConvertToDecimal(buffer, size, flags, digits,
352 rounding, Fortran::decimal::BinaryFloatingPointNumber<53>(x));
353 }
354
355 #if LDBL_MANT_DIG == 64
ConvertLongDoubleToDecimal(char * buffer,std::size_t size,enum DecimalConversionFlags flags,int digits,enum FortranRounding rounding,long double x)356 ConversionToDecimalResult ConvertLongDoubleToDecimal(char *buffer,
357 std::size_t size, enum DecimalConversionFlags flags, int digits,
358 enum FortranRounding rounding, long double x) {
359 return Fortran::decimal::ConvertToDecimal(buffer, size, flags, digits,
360 rounding, Fortran::decimal::BinaryFloatingPointNumber<64>(x));
361 }
362 #elif LDBL_MANT_DIG == 113
ConvertLongDoubleToDecimal(char * buffer,std::size_t size,enum DecimalConversionFlags flags,int digits,enum FortranRounding rounding,long double x)363 ConversionToDecimalResult ConvertLongDoubleToDecimal(char *buffer,
364 std::size_t size, enum DecimalConversionFlags flags, int digits,
365 enum FortranRounding rounding, long double x) {
366 return Fortran::decimal::ConvertToDecimal(buffer, size, flags, digits,
367 rounding, Fortran::decimal::BinaryFloatingPointNumber<113>(x));
368 }
369 #endif
370
371 RT_EXT_API_GROUP_END
372 } // extern "C"
373
374 template <int PREC, int LOG10RADIX>
375 template <typename STREAM>
Dump(STREAM & o) const376 STREAM &BigRadixFloatingPointNumber<PREC, LOG10RADIX>::Dump(STREAM &o) const {
377 if (isNegative_) {
378 o << '-';
379 }
380 o << "10**(" << exponent_ << ") * ... (rounding "
381 << static_cast<int>(rounding_) << ")\n";
382 for (int j{digits_}; --j >= 0;) {
383 std::string str{std::to_string(digit_[j])};
384 o << std::string(20 - str.size(), ' ') << str << " [" << j << ']';
385 if (j + 1 == digitLimit_) {
386 o << " (limit)";
387 }
388 o << '\n';
389 }
390 return o;
391 }
392 } // namespace Fortran::decimal
393