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>
BigRadixFloatingPointNumber(BinaryFloatingPointNumber<PREC> x,enum FortranRounding rounding)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
ConvertToDecimal(char * buffer,std::size_t n,enum DecimalConversionFlags flags,int maxDigits) const80 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 // Treat the MSD specially: don't emit leading zeroes.
104 Digit dig{digit_[digits_ - 1]};
105 char stack[LOG10RADIX], *sp{stack};
106 for (int k{0}; k < log10Radix; k += 2) {
107 Digit newDig{dig / 100};
108 auto d{static_cast<std::uint32_t>(dig) -
109 std::uint32_t{100} * static_cast<std::uint32_t>(newDig)};
110 dig = newDig;
111 const char *q{lut + d + d};
112 *sp++ = q[1];
113 *sp++ = q[0];
114 }
115 while (sp > stack && sp[-1] == '0') {
116 --sp;
117 }
118 while (sp > stack) {
119 *p++ = *--sp;
120 }
121 for (int j{digits_ - 1}; j-- > 0;) {
122 Digit dig{digit_[j]};
123 char *reverse{p += log10Radix};
124 for (int k{0}; k < log10Radix; k += 2) {
125 Digit newDig{dig / 100};
126 auto d{static_cast<std::uint32_t>(dig) -
127 std::uint32_t{100} * static_cast<std::uint32_t>(newDig)};
128 dig = newDig;
129 const char *q{lut + d + d};
130 *--reverse = q[1];
131 *--reverse = q[0];
132 }
133 }
134 // Adjust exponent so the effective decimal point is to
135 // the left of the first digit.
136 int expo = exponent_ + p - start;
137 // Trim trailing zeroes.
138 while (p[-1] == '0') {
139 --p;
140 }
141 char *end{start + maxDigits};
142 if (maxDigits == 0) {
143 p = end;
144 }
145 if (p <= end) {
146 *p = '\0';
147 return {buffer, static_cast<std::size_t>(p - buffer), expo, Exact};
148 } else {
149 // Apply a digit limit, possibly with rounding.
150 bool incr{false};
151 switch (rounding_) {
152 case RoundNearest:
153 incr = *end > '5' ||
154 (*end == '5' && (p > end + 1 || ((end[-1] - '0') & 1) != 0));
155 break;
156 case RoundUp:
157 incr = !isNegative_;
158 break;
159 case RoundDown:
160 incr = isNegative_;
161 break;
162 case RoundToZero:
163 break;
164 case RoundCompatible:
165 incr = *end >= '5';
166 break;
167 }
168 p = end;
169 if (incr) {
170 while (p > start && p[-1] == '9') {
171 --p;
172 }
173 if (p == start) {
174 *p++ = '1';
175 ++expo;
176 } else {
177 ++p[-1];
178 }
179 }
180
181 *p = '\0';
182 return {buffer, static_cast<std::size_t>(p - buffer), expo, Inexact};
183 }
184 }
185
186 template <int PREC, int LOG10RADIX>
Mean(const BigRadixFloatingPointNumber & that)187 bool BigRadixFloatingPointNumber<PREC, LOG10RADIX>::Mean(
188 const BigRadixFloatingPointNumber &that) {
189 while (digits_ < that.digits_) {
190 digit_[digits_++] = 0;
191 }
192 int carry{0};
193 for (int j{0}; j < that.digits_; ++j) {
194 Digit v{digit_[j] + that.digit_[j] + carry};
195 if (v >= radix) {
196 digit_[j] = v - radix;
197 carry = 1;
198 } else {
199 digit_[j] = v;
200 carry = 0;
201 }
202 }
203 if (carry != 0) {
204 AddCarry(that.digits_, carry);
205 }
206 return DivideBy<2>() != 0;
207 }
208
209 template <int PREC, int LOG10RADIX>
Minimize(BigRadixFloatingPointNumber && less,BigRadixFloatingPointNumber && more)210 void BigRadixFloatingPointNumber<PREC, LOG10RADIX>::Minimize(
211 BigRadixFloatingPointNumber &&less, BigRadixFloatingPointNumber &&more) {
212 int leastExponent{exponent_};
213 if (less.exponent_ < leastExponent) {
214 leastExponent = less.exponent_;
215 }
216 if (more.exponent_ < leastExponent) {
217 leastExponent = more.exponent_;
218 }
219 while (exponent_ > leastExponent) {
220 --exponent_;
221 MultiplyBy<10>();
222 }
223 while (less.exponent_ > leastExponent) {
224 --less.exponent_;
225 less.MultiplyBy<10>();
226 }
227 while (more.exponent_ > leastExponent) {
228 --more.exponent_;
229 more.MultiplyBy<10>();
230 }
231 if (less.Mean(*this)) {
232 less.AddCarry(); // round up
233 }
234 if (!more.Mean(*this)) {
235 more.Decrement(); // round down
236 }
237 while (less.digits_ < more.digits_) {
238 less.digit_[less.digits_++] = 0;
239 }
240 while (more.digits_ < less.digits_) {
241 more.digit_[more.digits_++] = 0;
242 }
243 int digits{more.digits_};
244 int same{0};
245 while (same < digits &&
246 less.digit_[digits - 1 - same] == more.digit_[digits - 1 - same]) {
247 ++same;
248 }
249 if (same == digits) {
250 return;
251 }
252 digits_ = same + 1;
253 int offset{digits - digits_};
254 exponent_ += offset * log10Radix;
255 for (int j{0}; j < digits_; ++j) {
256 digit_[j] = more.digit_[j + offset];
257 }
258 Digit least{less.digit_[offset]};
259 Digit my{digit_[0]};
260 while (true) {
261 Digit q{my / 10u};
262 Digit r{my - 10 * q};
263 Digit lq{least / 10u};
264 Digit lr{least - 10 * lq};
265 if (r != 0 && lq == q) {
266 Digit sub{(r - lr) >> 1};
267 digit_[0] -= sub;
268 break;
269 } else {
270 least = lq;
271 my = q;
272 DivideBy<10>();
273 ++exponent_;
274 }
275 }
276 Normalize();
277 }
278
279 template <int PREC>
ConvertToDecimal(char * buffer,std::size_t size,enum DecimalConversionFlags flags,int digits,enum FortranRounding rounding,BinaryFloatingPointNumber<PREC> x)280 ConversionToDecimalResult ConvertToDecimal(char *buffer, std::size_t size,
281 enum DecimalConversionFlags flags, int digits,
282 enum FortranRounding rounding, BinaryFloatingPointNumber<PREC> x) {
283 if (x.IsNaN()) {
284 return {"NaN", 3, 0, Invalid};
285 } else if (x.IsInfinite()) {
286 if (x.IsNegative()) {
287 return {"-Inf", 4, 0, Exact};
288 } else if (flags & AlwaysSign) {
289 return {"+Inf", 4, 0, Exact};
290 } else {
291 return {"Inf", 3, 0, Exact};
292 }
293 } else {
294 using Big = BigRadixFloatingPointNumber<PREC>;
295 Big number{x, rounding};
296 if ((flags & Minimize) && !x.IsZero()) {
297 // To emit the fewest decimal digits necessary to represent the value
298 // in such a way that decimal-to-binary conversion to the same format
299 // with a fixed assumption about rounding will return the same binary
300 // value, we also perform binary-to-decimal conversion on the two
301 // binary values immediately adjacent to this one, use them to identify
302 // the bounds of the range of decimal values that will map back to the
303 // original binary value, and find a (not necessary unique) shortest
304 // decimal sequence in that range.
305 using Binary = typename Big::Real;
306 Binary less{x};
307 less.Previous();
308 Binary more{x};
309 if (!x.IsMaximalFiniteMagnitude()) {
310 more.Next();
311 }
312 number.Minimize(Big{less, rounding}, Big{more, rounding});
313 } else {
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" {
ConvertFloatToDecimal(char * buffer,std::size_t size,enum DecimalConversionFlags flags,int digits,enum FortranRounding rounding,float x)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
ConvertDoubleToDecimal(char * buffer,std::size_t size,enum DecimalConversionFlags flags,int digits,enum FortranRounding rounding,double x)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 __x86_64__ && !defined(_MSC_VER)
ConvertLongDoubleToDecimal(char * buffer,std::size_t size,enum DecimalConversionFlags flags,int digits,enum FortranRounding rounding,long double x)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 #endif
361 }
362
363 template <int PREC, int LOG10RADIX>
364 template <typename STREAM>
Dump(STREAM & o) const365 STREAM &BigRadixFloatingPointNumber<PREC, LOG10RADIX>::Dump(STREAM &o) const {
366 if (isNegative_) {
367 o << '-';
368 }
369 o << "10**(" << exponent_ << ") * ...\n";
370 for (int j{digits_}; --j >= 0;) {
371 std::string str{std::to_string(digit_[j])};
372 o << std::string(20 - str.size(), ' ') << str << " [" << j << ']';
373 if (j + 1 == digitLimit_) {
374 o << " (limit)";
375 }
376 o << '\n';
377 }
378 return o;
379 }
380 } // namespace Fortran::decimal
381