NumeRe v1.1.4
NumeRe: Framework für Numerische Rechnungen
simple_decimal_conversion.h
Go to the documentation of this file.
1#ifndef FASTFLOAT_GENERIC_DECIMAL_TO_BINARY_H
2#define FASTFLOAT_GENERIC_DECIMAL_TO_BINARY_H
3
16#include "ascii_number.h"
17#include "decimal_to_binary.h"
18#include <cstdint>
19
20namespace fast_float {
21
22namespace detail {
23
24// remove all final zeroes
25inline void trim(decimal &h) {
26 while ((h.num_digits > 0) && (h.digits[h.num_digits - 1] == 0)) {
27 h.num_digits--;
28 }
29}
30
31
32
33inline uint32_t number_of_digits_decimal_left_shift(const decimal &h, uint32_t shift) {
34 shift &= 63;
35 constexpr uint16_t number_of_digits_decimal_left_shift_table[65] = {
36 0x0000, 0x0800, 0x0801, 0x0803, 0x1006, 0x1009, 0x100D, 0x1812, 0x1817,
37 0x181D, 0x2024, 0x202B, 0x2033, 0x203C, 0x2846, 0x2850, 0x285B, 0x3067,
38 0x3073, 0x3080, 0x388E, 0x389C, 0x38AB, 0x38BB, 0x40CC, 0x40DD, 0x40EF,
39 0x4902, 0x4915, 0x4929, 0x513E, 0x5153, 0x5169, 0x5180, 0x5998, 0x59B0,
40 0x59C9, 0x61E3, 0x61FD, 0x6218, 0x6A34, 0x6A50, 0x6A6D, 0x6A8B, 0x72AA,
41 0x72C9, 0x72E9, 0x7B0A, 0x7B2B, 0x7B4D, 0x8370, 0x8393, 0x83B7, 0x83DC,
42 0x8C02, 0x8C28, 0x8C4F, 0x9477, 0x949F, 0x94C8, 0x9CF2, 0x051C, 0x051C,
43 0x051C, 0x051C,
44 };
45 uint32_t x_a = number_of_digits_decimal_left_shift_table[shift];
46 uint32_t x_b = number_of_digits_decimal_left_shift_table[shift + 1];
47 uint32_t num_new_digits = x_a >> 11;
48 uint32_t pow5_a = 0x7FF & x_a;
49 uint32_t pow5_b = 0x7FF & x_b;
50 constexpr uint8_t
51 number_of_digits_decimal_left_shift_table_powers_of_5[0x051C] = {
52 5, 2, 5, 1, 2, 5, 6, 2, 5, 3, 1, 2, 5, 1, 5, 6, 2, 5, 7, 8, 1, 2, 5, 3,
53 9, 0, 6, 2, 5, 1, 9, 5, 3, 1, 2, 5, 9, 7, 6, 5, 6, 2, 5, 4, 8, 8, 2, 8,
54 1, 2, 5, 2, 4, 4, 1, 4, 0, 6, 2, 5, 1, 2, 2, 0, 7, 0, 3, 1, 2, 5, 6, 1,
55 0, 3, 5, 1, 5, 6, 2, 5, 3, 0, 5, 1, 7, 5, 7, 8, 1, 2, 5, 1, 5, 2, 5, 8,
56 7, 8, 9, 0, 6, 2, 5, 7, 6, 2, 9, 3, 9, 4, 5, 3, 1, 2, 5, 3, 8, 1, 4, 6,
57 9, 7, 2, 6, 5, 6, 2, 5, 1, 9, 0, 7, 3, 4, 8, 6, 3, 2, 8, 1, 2, 5, 9, 5,
58 3, 6, 7, 4, 3, 1, 6, 4, 0, 6, 2, 5, 4, 7, 6, 8, 3, 7, 1, 5, 8, 2, 0, 3,
59 1, 2, 5, 2, 3, 8, 4, 1, 8, 5, 7, 9, 1, 0, 1, 5, 6, 2, 5, 1, 1, 9, 2, 0,
60 9, 2, 8, 9, 5, 5, 0, 7, 8, 1, 2, 5, 5, 9, 6, 0, 4, 6, 4, 4, 7, 7, 5, 3,
61 9, 0, 6, 2, 5, 2, 9, 8, 0, 2, 3, 2, 2, 3, 8, 7, 6, 9, 5, 3, 1, 2, 5, 1,
62 4, 9, 0, 1, 1, 6, 1, 1, 9, 3, 8, 4, 7, 6, 5, 6, 2, 5, 7, 4, 5, 0, 5, 8,
63 0, 5, 9, 6, 9, 2, 3, 8, 2, 8, 1, 2, 5, 3, 7, 2, 5, 2, 9, 0, 2, 9, 8, 4,
64 6, 1, 9, 1, 4, 0, 6, 2, 5, 1, 8, 6, 2, 6, 4, 5, 1, 4, 9, 2, 3, 0, 9, 5,
65 7, 0, 3, 1, 2, 5, 9, 3, 1, 3, 2, 2, 5, 7, 4, 6, 1, 5, 4, 7, 8, 5, 1, 5,
66 6, 2, 5, 4, 6, 5, 6, 6, 1, 2, 8, 7, 3, 0, 7, 7, 3, 9, 2, 5, 7, 8, 1, 2,
67 5, 2, 3, 2, 8, 3, 0, 6, 4, 3, 6, 5, 3, 8, 6, 9, 6, 2, 8, 9, 0, 6, 2, 5,
68 1, 1, 6, 4, 1, 5, 3, 2, 1, 8, 2, 6, 9, 3, 4, 8, 1, 4, 4, 5, 3, 1, 2, 5,
69 5, 8, 2, 0, 7, 6, 6, 0, 9, 1, 3, 4, 6, 7, 4, 0, 7, 2, 2, 6, 5, 6, 2, 5,
70 2, 9, 1, 0, 3, 8, 3, 0, 4, 5, 6, 7, 3, 3, 7, 0, 3, 6, 1, 3, 2, 8, 1, 2,
71 5, 1, 4, 5, 5, 1, 9, 1, 5, 2, 2, 8, 3, 6, 6, 8, 5, 1, 8, 0, 6, 6, 4, 0,
72 6, 2, 5, 7, 2, 7, 5, 9, 5, 7, 6, 1, 4, 1, 8, 3, 4, 2, 5, 9, 0, 3, 3, 2,
73 0, 3, 1, 2, 5, 3, 6, 3, 7, 9, 7, 8, 8, 0, 7, 0, 9, 1, 7, 1, 2, 9, 5, 1,
74 6, 6, 0, 1, 5, 6, 2, 5, 1, 8, 1, 8, 9, 8, 9, 4, 0, 3, 5, 4, 5, 8, 5, 6,
75 4, 7, 5, 8, 3, 0, 0, 7, 8, 1, 2, 5, 9, 0, 9, 4, 9, 4, 7, 0, 1, 7, 7, 2,
76 9, 2, 8, 2, 3, 7, 9, 1, 5, 0, 3, 9, 0, 6, 2, 5, 4, 5, 4, 7, 4, 7, 3, 5,
77 0, 8, 8, 6, 4, 6, 4, 1, 1, 8, 9, 5, 7, 5, 1, 9, 5, 3, 1, 2, 5, 2, 2, 7,
78 3, 7, 3, 6, 7, 5, 4, 4, 3, 2, 3, 2, 0, 5, 9, 4, 7, 8, 7, 5, 9, 7, 6, 5,
79 6, 2, 5, 1, 1, 3, 6, 8, 6, 8, 3, 7, 7, 2, 1, 6, 1, 6, 0, 2, 9, 7, 3, 9,
80 3, 7, 9, 8, 8, 2, 8, 1, 2, 5, 5, 6, 8, 4, 3, 4, 1, 8, 8, 6, 0, 8, 0, 8,
81 0, 1, 4, 8, 6, 9, 6, 8, 9, 9, 4, 1, 4, 0, 6, 2, 5, 2, 8, 4, 2, 1, 7, 0,
82 9, 4, 3, 0, 4, 0, 4, 0, 0, 7, 4, 3, 4, 8, 4, 4, 9, 7, 0, 7, 0, 3, 1, 2,
83 5, 1, 4, 2, 1, 0, 8, 5, 4, 7, 1, 5, 2, 0, 2, 0, 0, 3, 7, 1, 7, 4, 2, 2,
84 4, 8, 5, 3, 5, 1, 5, 6, 2, 5, 7, 1, 0, 5, 4, 2, 7, 3, 5, 7, 6, 0, 1, 0,
85 0, 1, 8, 5, 8, 7, 1, 1, 2, 4, 2, 6, 7, 5, 7, 8, 1, 2, 5, 3, 5, 5, 2, 7,
86 1, 3, 6, 7, 8, 8, 0, 0, 5, 0, 0, 9, 2, 9, 3, 5, 5, 6, 2, 1, 3, 3, 7, 8,
87 9, 0, 6, 2, 5, 1, 7, 7, 6, 3, 5, 6, 8, 3, 9, 4, 0, 0, 2, 5, 0, 4, 6, 4,
88 6, 7, 7, 8, 1, 0, 6, 6, 8, 9, 4, 5, 3, 1, 2, 5, 8, 8, 8, 1, 7, 8, 4, 1,
89 9, 7, 0, 0, 1, 2, 5, 2, 3, 2, 3, 3, 8, 9, 0, 5, 3, 3, 4, 4, 7, 2, 6, 5,
90 6, 2, 5, 4, 4, 4, 0, 8, 9, 2, 0, 9, 8, 5, 0, 0, 6, 2, 6, 1, 6, 1, 6, 9,
91 4, 5, 2, 6, 6, 7, 2, 3, 6, 3, 2, 8, 1, 2, 5, 2, 2, 2, 0, 4, 4, 6, 0, 4,
92 9, 2, 5, 0, 3, 1, 3, 0, 8, 0, 8, 4, 7, 2, 6, 3, 3, 3, 6, 1, 8, 1, 6, 4,
93 0, 6, 2, 5, 1, 1, 1, 0, 2, 2, 3, 0, 2, 4, 6, 2, 5, 1, 5, 6, 5, 4, 0, 4,
94 2, 3, 6, 3, 1, 6, 6, 8, 0, 9, 0, 8, 2, 0, 3, 1, 2, 5, 5, 5, 5, 1, 1, 1,
95 5, 1, 2, 3, 1, 2, 5, 7, 8, 2, 7, 0, 2, 1, 1, 8, 1, 5, 8, 3, 4, 0, 4, 5,
96 4, 1, 0, 1, 5, 6, 2, 5, 2, 7, 7, 5, 5, 5, 7, 5, 6, 1, 5, 6, 2, 8, 9, 1,
97 3, 5, 1, 0, 5, 9, 0, 7, 9, 1, 7, 0, 2, 2, 7, 0, 5, 0, 7, 8, 1, 2, 5, 1,
98 3, 8, 7, 7, 7, 8, 7, 8, 0, 7, 8, 1, 4, 4, 5, 6, 7, 5, 5, 2, 9, 5, 3, 9,
99 5, 8, 5, 1, 1, 3, 5, 2, 5, 3, 9, 0, 6, 2, 5, 6, 9, 3, 8, 8, 9, 3, 9, 0,
100 3, 9, 0, 7, 2, 2, 8, 3, 7, 7, 6, 4, 7, 6, 9, 7, 9, 2, 5, 5, 6, 7, 6, 2,
101 6, 9, 5, 3, 1, 2, 5, 3, 4, 6, 9, 4, 4, 6, 9, 5, 1, 9, 5, 3, 6, 1, 4, 1,
102 8, 8, 8, 2, 3, 8, 4, 8, 9, 6, 2, 7, 8, 3, 8, 1, 3, 4, 7, 6, 5, 6, 2, 5,
103 1, 7, 3, 4, 7, 2, 3, 4, 7, 5, 9, 7, 6, 8, 0, 7, 0, 9, 4, 4, 1, 1, 9, 2,
104 4, 4, 8, 1, 3, 9, 1, 9, 0, 6, 7, 3, 8, 2, 8, 1, 2, 5, 8, 6, 7, 3, 6, 1,
105 7, 3, 7, 9, 8, 8, 4, 0, 3, 5, 4, 7, 2, 0, 5, 9, 6, 2, 2, 4, 0, 6, 9, 5,
106 9, 5, 3, 3, 6, 9, 1, 4, 0, 6, 2, 5,
107 };
108 const uint8_t *pow5 =
109 &number_of_digits_decimal_left_shift_table_powers_of_5[pow5_a];
110 uint32_t i = 0;
111 uint32_t n = pow5_b - pow5_a;
112 for (; i < n; i++) {
113 if (i >= h.num_digits) {
114 return num_new_digits - 1;
115 } else if (h.digits[i] == pow5[i]) {
116 continue;
117 } else if (h.digits[i] < pow5[i]) {
118 return num_new_digits - 1;
119 } else {
120 return num_new_digits;
121 }
122 }
123 return num_new_digits;
124}
125
126inline uint64_t round(decimal &h) {
127 if ((h.num_digits == 0) || (h.decimal_point < 0)) {
128 return 0;
129 } else if (h.decimal_point > 18) {
130 return UINT64_MAX;
131 }
132 // at this point, we know that h.decimal_point >= 0
133 uint32_t dp = uint32_t(h.decimal_point);
134 uint64_t n = 0;
135 for (uint32_t i = 0; i < dp; i++) {
136 n = (10 * n) + ((i < h.num_digits) ? h.digits[i] : 0);
137 }
138 bool round_up = false;
139 if (dp < h.num_digits) {
140 round_up = h.digits[dp] >= 5; // normally, we round up
141 // but we may need to round to even!
142 if ((h.digits[dp] == 5) && (dp + 1 == h.num_digits)) {
143 round_up = h.truncated || ((dp > 0) && (1 & h.digits[dp - 1]));
144 }
145 }
146 if (round_up) {
147 n++;
148 }
149 return n;
150}
151
152// computes h * 2^-shift
153inline void decimal_left_shift(decimal &h, uint32_t shift) {
154 if (h.num_digits == 0) {
155 return;
156 }
157 uint32_t num_new_digits = number_of_digits_decimal_left_shift(h, shift);
158 int32_t read_index = int32_t(h.num_digits - 1);
159 uint32_t write_index = h.num_digits - 1 + num_new_digits;
160 uint64_t n = 0;
161
162 while (read_index >= 0) {
163 n += uint64_t(h.digits[read_index]) << shift;
164 uint64_t quotient = n / 10;
165 uint64_t remainder = n - (10 * quotient);
166 if (write_index < max_digits) {
167 h.digits[write_index] = uint8_t(remainder);
168 } else if (remainder > 0) {
169 h.truncated = true;
170 }
171 n = quotient;
172 write_index--;
173 read_index--;
174 }
175 while (n > 0) {
176 uint64_t quotient = n / 10;
177 uint64_t remainder = n - (10 * quotient);
178 if (write_index < max_digits) {
179 h.digits[write_index] = uint8_t(remainder);
180 } else if (remainder > 0) {
181 h.truncated = true;
182 }
183 n = quotient;
184 write_index--;
185 }
186 h.num_digits += num_new_digits;
187 if (h.num_digits > max_digits) {
188 h.num_digits = max_digits;
189 }
190 h.decimal_point += int32_t(num_new_digits);
191 trim(h);
192}
193
194// computes h * 2^shift
195inline void decimal_right_shift(decimal &h, uint32_t shift) {
196 uint32_t read_index = 0;
197 uint32_t write_index = 0;
198
199 uint64_t n = 0;
200
201 while ((n >> shift) == 0) {
202 if (read_index < h.num_digits) {
203 n = (10 * n) + h.digits[read_index++];
204 } else if (n == 0) {
205 return;
206 } else {
207 while ((n >> shift) == 0) {
208 n = 10 * n;
209 read_index++;
210 }
211 break;
212 }
213 }
214 h.decimal_point -= int32_t(read_index - 1);
215 if (h.decimal_point < -decimal_point_range) { // it is zero
216 h.num_digits = 0;
217 h.decimal_point = 0;
218 h.negative = false;
219 h.truncated = false;
220 return;
221 }
222 uint64_t mask = (uint64_t(1) << shift) - 1;
223 while (read_index < h.num_digits) {
224 uint8_t new_digit = uint8_t(n >> shift);
225 n = (10 * (n & mask)) + h.digits[read_index++];
226 h.digits[write_index++] = new_digit;
227 }
228 while (n > 0) {
229 uint8_t new_digit = uint8_t(n >> shift);
230 n = 10 * (n & mask);
231 if (write_index < max_digits) {
232 h.digits[write_index++] = new_digit;
233 } else if (new_digit > 0) {
234 h.truncated = true;
235 }
236 }
237 h.num_digits = write_index;
238 trim(h);
239}
240
241} // namespace detail
242
243template <typename binary>
245 adjusted_mantissa answer;
246 if (d.num_digits == 0) {
247 // should be zero
248 answer.power2 = 0;
249 answer.mantissa = 0;
250 return answer;
251 }
252 // At this point, going further, we can assume that d.num_digits > 0.
253 //
254 // We want to guard against excessive decimal point values because
255 // they can result in long running times. Indeed, we do
256 // shifts by at most 60 bits. We have that log(10**400)/log(2**60) ~= 22
257 // which is fine, but log(10**299995)/log(2**60) ~= 16609 which is not
258 // fine (runs for a long time).
259 //
260 if(d.decimal_point < -324) {
261 // We have something smaller than 1e-324 which is always zero
262 // in binary64 and binary32.
263 // It should be zero.
264 answer.power2 = 0;
265 answer.mantissa = 0;
266 return answer;
267 } else if(d.decimal_point >= 310) {
268 // We have something at least as large as 0.1e310 which is
269 // always infinite.
270 answer.power2 = binary::infinite_power();
271 answer.mantissa = 0;
272 return answer;
273 }
274 constexpr uint32_t max_shift = 60;
275 constexpr uint32_t num_powers = 19;
276 constexpr uint8_t decimal_powers[19] = {
277 0, 3, 6, 9, 13, 16, 19, 23, 26, 29, //
278 33, 36, 39, 43, 46, 49, 53, 56, 59, //
279 };
280 int32_t exp2 = 0;
281 while (d.decimal_point > 0) {
282 uint32_t n = uint32_t(d.decimal_point);
283 uint32_t shift = (n < num_powers) ? decimal_powers[n] : max_shift;
285 if (d.decimal_point < -decimal_point_range) {
286 // should be zero
287 answer.power2 = 0;
288 answer.mantissa = 0;
289 return answer;
290 }
291 exp2 += int32_t(shift);
292 }
293 // We shift left toward [1/2 ... 1].
294 while (d.decimal_point <= 0) {
295 uint32_t shift;
296 if (d.decimal_point == 0) {
297 if (d.digits[0] >= 5) {
298 break;
299 }
300 shift = (d.digits[0] < 2) ? 2 : 1;
301 } else {
302 uint32_t n = uint32_t(-d.decimal_point);
303 shift = (n < num_powers) ? decimal_powers[n] : max_shift;
304 }
306 if (d.decimal_point > decimal_point_range) {
307 // we want to get infinity:
308 answer.power2 = binary::infinite_power();
309 answer.mantissa = 0;
310 return answer;
311 }
312 exp2 -= int32_t(shift);
313 }
314 // We are now in the range [1/2 ... 1] but the binary format uses [1 ... 2].
315 exp2--;
316 constexpr int32_t minimum_exponent = binary::minimum_exponent();
317 while ((minimum_exponent + 1) > exp2) {
318 uint32_t n = uint32_t((minimum_exponent + 1) - exp2);
319 if (n > max_shift) {
320 n = max_shift;
321 }
323 exp2 += int32_t(n);
324 }
325 if ((exp2 - minimum_exponent) >= binary::infinite_power()) {
326 answer.power2 = binary::infinite_power();
327 answer.mantissa = 0;
328 return answer;
329 }
330
331 const int mantissa_size_in_bits = binary::mantissa_explicit_bits() + 1;
332 detail::decimal_left_shift(d, mantissa_size_in_bits);
333
334 uint64_t mantissa = detail::round(d);
335 // It is possible that we have an overflow, in which case we need
336 // to shift back.
337 if(mantissa >= (uint64_t(1) << mantissa_size_in_bits)) {
339 exp2 += 1;
340 mantissa = detail::round(d);
341 if ((exp2 - minimum_exponent) >= binary::infinite_power()) {
342 answer.power2 = binary::infinite_power();
343 answer.mantissa = 0;
344 return answer;
345 }
346 }
347 answer.power2 = exp2 - binary::minimum_exponent();
348 if(mantissa < (uint64_t(1) << binary::mantissa_explicit_bits())) { answer.power2--; }
349 answer.mantissa = mantissa & ((uint64_t(1) << binary::mantissa_explicit_bits()) - 1);
350 return answer;
351}
352
353template <typename binary>
354adjusted_mantissa parse_long_mantissa(const char *first, const char* last, parse_options options) {
355 decimal d = parse_decimal(first, last, options);
356 return compute_float<binary>(d);
357}
358
359} // namespace fast_float
360#endif
CONSTDATA date::last_spec last
Definition: date.h:1989
uint32_t number_of_digits_decimal_left_shift(const decimal &h, uint32_t shift)
void decimal_left_shift(decimal &h, uint32_t shift)
uint64_t round(decimal &h)
void decimal_right_shift(decimal &h, uint32_t shift)
adjusted_mantissa parse_long_mantissa(const char *first, const char *last, parse_options options)
fastfloat_really_inline adjusted_mantissa compute_float(int64_t q, uint64_t w) noexcept