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 | RT_EXT_API_GROUP_BEGIN |
340 | |
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 | |
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 |
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 |
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> |
376 | 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 | |