1// Copyright 2020-2023 Daniel Lemire
2// Copyright 2023 Matt Borland
3// Distributed under the Boost Software License, Version 1.0.
4// https://www.boost.org/LICENSE_1_0.txt
5//
6// Derivative of: https://github.com/fastfloat/fast_float
7
8#ifndef BOOST_CHARCONV_DETAIL_FASTFLOAT_DIGIT_COMPARISON_HPP
9#define BOOST_CHARCONV_DETAIL_FASTFLOAT_DIGIT_COMPARISON_HPP
10
11#include <boost/charconv/detail/fast_float/float_common.hpp>
12#include <boost/charconv/detail/fast_float/bigint.hpp>
13#include <boost/charconv/detail/fast_float/ascii_number.hpp>
14#include <algorithm>
15#include <cstdint>
16#include <cstring>
17#include <iterator>
18
19namespace boost { namespace charconv { namespace detail { namespace fast_float {
20
21// 1e0 to 1e19
22constexpr static uint64_t powers_of_ten_uint64[] = {
23 1UL, 10UL, 100UL, 1000UL, 10000UL, 100000UL, 1000000UL, 10000000UL, 100000000UL,
24 1000000000UL, 10000000000UL, 100000000000UL, 1000000000000UL, 10000000000000UL,
25 100000000000000UL, 1000000000000000UL, 10000000000000000UL, 100000000000000000UL,
26 1000000000000000000UL, 10000000000000000000UL};
27
28// calculate the exponent, in scientific notation, of the number.
29// this algorithm is not even close to optimized, but it has no practical
30// effect on performance: in order to have a faster algorithm, we'd need
31// to slow down performance for faster algorithms, and this is still fast.
32template <typename UC>
33BOOST_FORCEINLINE BOOST_CHARCONV_FASTFLOAT_CONSTEXPR14
34int32_t scientific_exponent(parsed_number_string_t<UC> & num) noexcept {
35 uint64_t mantissa = num.mantissa;
36 int32_t exponent = int32_t(num.exponent);
37 while (mantissa >= 10000) {
38 mantissa /= 10000;
39 exponent += 4;
40 }
41 while (mantissa >= 100) {
42 mantissa /= 100;
43 exponent += 2;
44 }
45 while (mantissa >= 10) {
46 mantissa /= 10;
47 exponent += 1;
48 }
49 return exponent;
50}
51
52// this converts a native floating-point number to an extended-precision float.
53template <typename T>
54BOOST_FORCEINLINE BOOST_CHARCONV_FASTFLOAT_CONSTEXPR20
55adjusted_mantissa to_extended(T value) noexcept {
56 using equiv_uint = typename binary_format<T>::equiv_uint;
57 constexpr equiv_uint exponent_mask = binary_format<T>::exponent_mask();
58 constexpr equiv_uint mantissa_mask = binary_format<T>::mantissa_mask();
59 constexpr equiv_uint hidden_bit_mask = binary_format<T>::hidden_bit_mask();
60
61 adjusted_mantissa am;
62 int32_t bias = binary_format<T>::mantissa_explicit_bits() - binary_format<T>::minimum_exponent();
63 equiv_uint bits;
64#if BOOST_CHARCONV_FASTFLOAT_HAS_BIT_CAST
65 bits = std::bit_cast<equiv_uint>(value);
66#else
67 ::memcpy(dest: &bits, src: &value, n: sizeof(T));
68#endif
69 if ((bits & exponent_mask) == 0) {
70 // denormal
71 am.power2 = 1 - bias;
72 am.mantissa = bits & mantissa_mask;
73 } else {
74 // normal
75 am.power2 = int32_t((bits & exponent_mask) >> binary_format<T>::mantissa_explicit_bits());
76 am.power2 -= bias;
77 am.mantissa = (bits & mantissa_mask) | hidden_bit_mask;
78 }
79
80 return am;
81}
82
83// get the extended precision value of the halfway point between b and b+u.
84// we are given a native float that represents b, so we need to adjust it
85// halfway between b and b+u.
86template <typename T>
87BOOST_FORCEINLINE BOOST_CHARCONV_FASTFLOAT_CONSTEXPR20
88adjusted_mantissa to_extended_halfway(T value) noexcept {
89 adjusted_mantissa am = to_extended(value);
90 am.mantissa <<= 1;
91 am.mantissa += 1;
92 am.power2 -= 1;
93 return am;
94}
95
96// round an extended-precision float to the nearest machine float.
97template <typename T, typename callback>
98BOOST_FORCEINLINE BOOST_CHARCONV_FASTFLOAT_CONSTEXPR14
99void round(adjusted_mantissa& am, callback cb) noexcept {
100 int32_t mantissa_shift = 64 - binary_format<T>::mantissa_explicit_bits() - 1;
101 if (-am.power2 >= mantissa_shift) {
102 // have a denormal float
103 int32_t shift = -am.power2 + 1;
104 cb(am, std::min<int32_t>(a: shift, b: 64));
105 // check for round-up: if rounding-nearest carried us to the hidden bit.
106 am.power2 = (am.mantissa < (uint64_t(1) << binary_format<T>::mantissa_explicit_bits())) ? 0 : 1;
107 return;
108 }
109
110 // have a normal float, use the default shift.
111 cb(am, mantissa_shift);
112
113 // check for carry
114 if (am.mantissa >= (uint64_t(2) << binary_format<T>::mantissa_explicit_bits())) {
115 am.mantissa = (uint64_t(1) << binary_format<T>::mantissa_explicit_bits());
116 am.power2++;
117 }
118
119 // check for infinite: we could have carried to an infinite power
120 am.mantissa &= ~(uint64_t(1) << binary_format<T>::mantissa_explicit_bits());
121 if (am.power2 >= binary_format<T>::infinite_power()) {
122 am.power2 = binary_format<T>::infinite_power();
123 am.mantissa = 0;
124 }
125}
126
127template <typename callback>
128BOOST_FORCEINLINE BOOST_CHARCONV_FASTFLOAT_CONSTEXPR14
129void round_nearest_tie_even(adjusted_mantissa& am, int32_t shift, callback cb) noexcept {
130 const uint64_t mask
131 = (shift == 64)
132 ? UINT64_MAX
133 : (uint64_t(1) << shift) - 1;
134 const uint64_t halfway
135 = (shift == 0)
136 ? 0
137 : uint64_t(1) << (shift - 1);
138 uint64_t truncated_bits = am.mantissa & mask;
139 bool is_above = truncated_bits > halfway;
140 bool is_halfway = truncated_bits == halfway;
141
142 // shift digits into position
143 if (shift == 64) {
144 am.mantissa = 0;
145 } else {
146 am.mantissa >>= shift;
147 }
148 am.power2 += shift;
149
150 bool is_odd = (am.mantissa & 1) == 1;
151 am.mantissa += uint64_t(cb(is_odd, is_halfway, is_above));
152}
153
154BOOST_FORCEINLINE BOOST_CHARCONV_FASTFLOAT_CONSTEXPR14
155void round_down(adjusted_mantissa& am, int32_t shift) noexcept {
156 if (shift == 64) {
157 am.mantissa = 0;
158 } else {
159 am.mantissa >>= shift;
160 }
161 am.power2 += shift;
162}
163template <typename UC>
164BOOST_FORCEINLINE BOOST_CHARCONV_FASTFLOAT_CONSTEXPR20
165void skip_zeros(UC const * & first, UC const * last) noexcept {
166 uint64_t val;
167 while (!cpp20_and_in_constexpr() && std::distance(first, last) >= int_cmp_len<UC>()) {
168 ::memcpy(dest: &val, src: first, n: sizeof(uint64_t));
169 if (val != int_cmp_zeros<UC>()) {
170 break;
171 }
172 first += int_cmp_len<UC>();
173 }
174 while (first != last) {
175 if (*first != UC('0')) {
176 break;
177 }
178 first++;
179 }
180}
181
182// determine if any non-zero digits were truncated.
183// all characters must be valid digits.
184template <typename UC>
185BOOST_FORCEINLINE BOOST_CHARCONV_FASTFLOAT_CONSTEXPR20
186bool is_truncated(UC const * first, UC const * last) noexcept {
187 // do 8-bit optimizations, can just compare to 8 literal 0s.
188 uint64_t val;
189 while (!cpp20_and_in_constexpr() && std::distance(first, last) >= int_cmp_len<UC>()) {
190 ::memcpy(dest: &val, src: first, n: sizeof(uint64_t));
191 if (val != int_cmp_zeros<UC>()) {
192 return true;
193 }
194 first += int_cmp_len<UC>();
195 }
196 while (first != last) {
197 if (*first != UC('0')) {
198 return true;
199 }
200 ++first;
201 }
202 return false;
203}
204template <typename UC>
205BOOST_FORCEINLINE BOOST_CHARCONV_FASTFLOAT_CONSTEXPR20
206bool is_truncated(span<const UC> s) noexcept {
207 return is_truncated(s.ptr, s.ptr + s.len());
208}
209
210BOOST_FORCEINLINE BOOST_CHARCONV_FASTFLOAT_CONSTEXPR20
211void parse_eight_digits(const char16_t*& , limb& , size_t& , size_t& ) noexcept {
212 // currently unused
213}
214
215BOOST_FORCEINLINE BOOST_CHARCONV_FASTFLOAT_CONSTEXPR20
216void parse_eight_digits(const char32_t*& , limb& , size_t& , size_t& ) noexcept {
217 // currently unused
218}
219
220BOOST_FORCEINLINE BOOST_CHARCONV_FASTFLOAT_CONSTEXPR20
221void parse_eight_digits(const char*& p, limb& value, size_t& counter, size_t& count) noexcept {
222 value = value * 100000000 + parse_eight_digits_unrolled(chars: p);
223 p += 8;
224 counter += 8;
225 count += 8;
226}
227
228template <typename UC>
229BOOST_FORCEINLINE BOOST_CHARCONV_FASTFLOAT_CONSTEXPR14
230void parse_one_digit(UC const *& p, limb& value, size_t& counter, size_t& count) noexcept {
231 value = value * 10 + limb(*p - UC('0'));
232 p++;
233 counter++;
234 count++;
235}
236
237BOOST_FORCEINLINE BOOST_CHARCONV_FASTFLOAT_CONSTEXPR20
238void add_native(bigint& big, limb power, limb value) noexcept {
239 big.mul(y: power);
240 big.add(y: value);
241}
242
243BOOST_FORCEINLINE BOOST_CHARCONV_FASTFLOAT_CONSTEXPR20
244void round_up_bigint(bigint& big, size_t& count) noexcept {
245 // need to round-up the digits, but need to avoid rounding
246 // ....9999 to ...10000, which could cause a false halfway point.
247 add_native(big, power: 10, value: 1);
248 count++;
249}
250
251// parse the significant digits into a big integer
252template <typename UC>
253inline BOOST_CHARCONV_FASTFLOAT_CONSTEXPR20
254void parse_mantissa(bigint& result, parsed_number_string_t<UC>& num, size_t max_digits, size_t& digits) noexcept {
255 // try to minimize the number of big integer and scalar multiplication.
256 // therefore, try to parse 8 digits at a time, and multiply by the largest
257 // scalar value (9 or 19 digits) for each step.
258 size_t counter = 0;
259 digits = 0;
260 limb value = 0;
261#ifdef BOOST_CHARCONV_FASTFLOAT_64BIT_LIMB
262 constexpr size_t step = 19;
263#else
264 constexpr size_t step = 9;
265#endif
266
267 // process all integer digits.
268 UC const * p = num.integer.ptr;
269 UC const * pend = p + num.integer.len();
270 skip_zeros(p, pend);
271 // process all digits, in increments of step per loop
272 while (p != pend) {
273 if (std::is_same<UC,char>::value) {
274 while ((std::distance(p, pend) >= 8) && (step - counter >= 8) && (max_digits - digits >= 8)) {
275 parse_eight_digits(p, value, counter, digits);
276 }
277 }
278 while (counter < step && p != pend && digits < max_digits) {
279 parse_one_digit(p, value, counter, digits);
280 }
281 if (digits == max_digits) {
282 // add the temporary value, then check if we've truncated any digits
283 add_native(big&: result, power: limb(powers_of_ten_uint64[counter]), value);
284 bool truncated = is_truncated(p, pend);
285 if (num.fraction.ptr != nullptr) {
286 truncated |= is_truncated(num.fraction);
287 }
288 if (truncated) {
289 round_up_bigint(big&: result, count&: digits);
290 }
291 return;
292 } else {
293 add_native(big&: result, power: limb(powers_of_ten_uint64[counter]), value);
294 counter = 0;
295 value = 0;
296 }
297 }
298
299 // add our fraction digits, if they're available.
300 if (num.fraction.ptr != nullptr) {
301 p = num.fraction.ptr;
302 pend = p + num.fraction.len();
303 if (digits == 0) {
304 skip_zeros(p, pend);
305 }
306 // process all digits, in increments of step per loop
307 while (p != pend) {
308 if (std::is_same<UC,char>::value) {
309 while ((std::distance(p, pend) >= 8) && (step - counter >= 8) && (max_digits - digits >= 8)) {
310 parse_eight_digits(p, value, counter, digits);
311 }
312 }
313 while (counter < step && p != pend && digits < max_digits) {
314 parse_one_digit(p, value, counter, digits);
315 }
316 if (digits == max_digits) {
317 // add the temporary value, then check if we've truncated any digits
318 add_native(big&: result, power: limb(powers_of_ten_uint64[counter]), value);
319 bool truncated = is_truncated(p, pend);
320 if (truncated) {
321 round_up_bigint(big&: result, count&: digits);
322 }
323 return;
324 } else {
325 add_native(big&: result, power: limb(powers_of_ten_uint64[counter]), value);
326 counter = 0;
327 value = 0;
328 }
329 }
330 }
331
332 if (counter != 0) {
333 add_native(big&: result, power: limb(powers_of_ten_uint64[counter]), value);
334 }
335}
336
337template <typename T>
338inline BOOST_CHARCONV_FASTFLOAT_CONSTEXPR20
339adjusted_mantissa positive_digit_comp(bigint& bigmant, int32_t exponent) noexcept {
340 BOOST_CHARCONV_FASTFLOAT_ASSERT(bigmant.pow10(uint32_t(exponent)));
341 adjusted_mantissa answer;
342 bool truncated;
343 answer.mantissa = bigmant.hi64(truncated);
344 int bias = binary_format<T>::mantissa_explicit_bits() - binary_format<T>::minimum_exponent();
345 answer.power2 = bigmant.bit_length() - 64 + bias;
346
347 round<T>(answer, [truncated](adjusted_mantissa& a, int32_t shift) {
348 round_nearest_tie_even(a, shift, [truncated](bool is_odd, bool is_halfway, bool is_above) -> bool {
349 return is_above || (is_halfway && truncated) || (is_odd && is_halfway);
350 });
351 });
352
353 return answer;
354}
355
356// the scaling here is quite simple: we have, for the real digits `m * 10^e`,
357// and for the theoretical digits `n * 2^f`. Since `e` is always negative,
358// to scale them identically, we do `n * 2^f * 5^-f`, so we now have `m * 2^e`.
359// we then need to scale by `2^(f- e)`, and then the two significant digits
360// are of the same magnitude.
361template <typename T>
362inline BOOST_CHARCONV_FASTFLOAT_CONSTEXPR20
363adjusted_mantissa negative_digit_comp(bigint& bigmant, adjusted_mantissa am, int32_t exponent) noexcept {
364 bigint& real_digits = bigmant;
365 int32_t real_exp = exponent;
366
367 // get the value of `b`, rounded down, and get a bigint representation of b+h
368 adjusted_mantissa am_b = am;
369 // gcc7 buf: use a lambda to remove the noexcept qualifier bug with -Wnoexcept-type.
370 round<T>(am_b, [](adjusted_mantissa&a, int32_t shift) { round_down(am&: a, shift); });
371 T b;
372 to_float(false, am_b, b);
373 adjusted_mantissa theor = to_extended_halfway(b);
374 bigint theor_digits(theor.mantissa);
375 int32_t theor_exp = theor.power2;
376
377 // scale real digits and theor digits to be same power.
378 int32_t pow2_exp = theor_exp - real_exp;
379 uint32_t pow5_exp = uint32_t(-real_exp);
380 if (pow5_exp != 0) {
381 BOOST_CHARCONV_FASTFLOAT_ASSERT(theor_digits.pow5(pow5_exp));
382 }
383 if (pow2_exp > 0) {
384 BOOST_CHARCONV_FASTFLOAT_ASSERT(theor_digits.pow2(uint32_t(pow2_exp)));
385 } else if (pow2_exp < 0) {
386 BOOST_CHARCONV_FASTFLOAT_ASSERT(real_digits.pow2(uint32_t(-pow2_exp)));
387 }
388
389 // compare digits, and use it to director rounding
390 int ord = real_digits.compare(other: theor_digits);
391 adjusted_mantissa answer = am;
392 round<T>(answer, [ord](adjusted_mantissa& a, int32_t shift) {
393 round_nearest_tie_even(a, shift, [ord](bool is_odd, bool, bool) -> bool {
394 if (ord > 0) {
395 return true;
396 } else if (ord < 0) {
397 return false;
398 } else {
399 return is_odd;
400 }
401 });
402 });
403
404 return answer;
405}
406
407// parse the significant digits as a big integer to unambiguously round
408// the significant digits. here, we are trying to determine how to round
409// an extended float representation close to `b+h`, halfway between `b`
410// (the float rounded-down) and `b+u`, the next positive float. this
411// algorithm is always correct, and uses one of two approaches. when
412// the exponent is positive relative to the significant digits (such as
413// 1234), we create a big-integer representation, get the high 64-bits,
414// determine if any lower bits are truncated, and use that to direct
415// rounding. in case of a negative exponent relative to the significant
416// digits (such as 1.2345), we create a theoretical representation of
417// `b` as a big-integer type, scaled to the same binary exponent as
418// the actual digits. we then compare the big integer representations
419// of both, and use that to direct rounding.
420template <typename T, typename UC>
421inline BOOST_CHARCONV_FASTFLOAT_CONSTEXPR20
422adjusted_mantissa digit_comp(parsed_number_string_t<UC>& num, adjusted_mantissa am) noexcept {
423 // remove the invalid exponent bias
424 am.power2 -= invalid_am_bias;
425
426 int32_t sci_exp = scientific_exponent(num);
427 size_t max_digits = binary_format<T>::max_digits();
428 size_t digits = 0;
429 bigint bigmant;
430 parse_mantissa(bigmant, num, max_digits, digits);
431 // can't underflow, since digits is at most max_digits.
432 int32_t exponent = sci_exp + 1 - int32_t(digits);
433 if (exponent >= 0) {
434 return positive_digit_comp<T>(bigmant, exponent);
435 } else {
436 return negative_digit_comp<T>(bigmant, am, exponent);
437 }
438}
439
440}}}} // namespace fast_float
441
442#endif
443

source code of boost/libs/charconv/include/boost/charconv/detail/fast_float/digit_comparison.hpp