94.74% Lines (54/57) 0.00% Functions (0/0)
TLA Baseline Branch
Line Hits Code Line Hits Code
1   // Copyright 2020-2023 Daniel Lemire 1   // Copyright 2020-2023 Daniel Lemire
2   // Copyright 2023 Matt Borland 2   // Copyright 2023 Matt Borland
3   // Distributed under the Boost Software License, Version 1.0. 3   // Distributed under the Boost Software License, Version 1.0.
4   // https://www.boost.org/LICENSE_1_0.txt 4   // https://www.boost.org/LICENSE_1_0.txt
5   // 5   //
6   // Derivative of: https://github.com/fastfloat/fast_float 6   // Derivative of: https://github.com/fastfloat/fast_float
7   7  
8   #ifndef BOOST_JSON_DETAIL_CHARCONV_DETAIL_FASTFLOAT_DECIMAL_TO_BINARY_HPP 8   #ifndef BOOST_JSON_DETAIL_CHARCONV_DETAIL_FASTFLOAT_DECIMAL_TO_BINARY_HPP
9   #define BOOST_JSON_DETAIL_CHARCONV_DETAIL_FASTFLOAT_DECIMAL_TO_BINARY_HPP 9   #define BOOST_JSON_DETAIL_CHARCONV_DETAIL_FASTFLOAT_DECIMAL_TO_BINARY_HPP
10   10  
11   #include <boost/json/detail/charconv/detail/fast_float/float_common.hpp> 11   #include <boost/json/detail/charconv/detail/fast_float/float_common.hpp>
12   #include <boost/json/detail/charconv/detail/fast_float/fast_table.hpp> 12   #include <boost/json/detail/charconv/detail/fast_float/fast_table.hpp>
13   #include <cfloat> 13   #include <cfloat>
14   #include <cinttypes> 14   #include <cinttypes>
15   #include <cmath> 15   #include <cmath>
16   #include <cstdint> 16   #include <cstdint>
17   #include <cstdlib> 17   #include <cstdlib>
18   #include <cstring> 18   #include <cstring>
19   19  
20   namespace boost { namespace json { namespace detail { namespace charconv { namespace detail { namespace fast_float { 20   namespace boost { namespace json { namespace detail { namespace charconv { namespace detail { namespace fast_float {
21   21  
22   // This will compute or rather approximate w * 5**q and return a pair of 64-bit words approximating 22   // This will compute or rather approximate w * 5**q and return a pair of 64-bit words approximating
23   // the result, with the "high" part corresponding to the most significant bits and the 23   // the result, with the "high" part corresponding to the most significant bits and the
24   // low part corresponding to the least significant bits. 24   // low part corresponding to the least significant bits.
25   // 25   //
26   template <int bit_precision> 26   template <int bit_precision>
27   BOOST_FORCEINLINE BOOST_JSON_FASTFLOAT_CONSTEXPR20 27   BOOST_FORCEINLINE BOOST_JSON_FASTFLOAT_CONSTEXPR20
28   value128 compute_product_approximation(int64_t q, uint64_t w) { 28   value128 compute_product_approximation(int64_t q, uint64_t w) {
HITCBC 29   2005879 const int index = 2 * int(q - powers::smallest_power_of_five); 29   2005879 const int index = 2 * int(q - powers::smallest_power_of_five);
30   // For small values of q, e.g., q in [0,27], the answer is always exact because 30   // For small values of q, e.g., q in [0,27], the answer is always exact because
31   // The line value128 firstproduct = full_multiplication(w, power_of_five_128[index]); 31   // The line value128 firstproduct = full_multiplication(w, power_of_five_128[index]);
32   // gives the exact answer. 32   // gives the exact answer.
HITCBC 33   4011758 value128 firstproduct = full_multiplication(w, powers::power_of_five_128[index]); 33   4011758 value128 firstproduct = full_multiplication(w, powers::power_of_five_128[index]);
34   static_assert((bit_precision >= 0) && (bit_precision <= 64), " precision should be in (0,64]"); 34   static_assert((bit_precision >= 0) && (bit_precision <= 64), " precision should be in (0,64]");
HITCBC 35   2005879 constexpr uint64_t precision_mask = (bit_precision < 64) ? 35   2005879 constexpr uint64_t precision_mask = (bit_precision < 64) ?
36   (uint64_t(0xFFFFFFFFFFFFFFFF) >> bit_precision) 36   (uint64_t(0xFFFFFFFFFFFFFFFF) >> bit_precision)
37   : uint64_t(0xFFFFFFFFFFFFFFFF); 37   : uint64_t(0xFFFFFFFFFFFFFFFF);
HITCBC 38   2005879 if((firstproduct.high & precision_mask) == precision_mask) { // could further guard with (lower + w < lower) 38   2005879 if((firstproduct.high & precision_mask) == precision_mask) { // could further guard with (lower + w < lower)
39   // regarding the second product, we only need secondproduct.high, but our expectation is that the compiler will optimize this extra work away if needed. 39   // regarding the second product, we only need secondproduct.high, but our expectation is that the compiler will optimize this extra work away if needed.
HITCBC 40   4995 value128 secondproduct = full_multiplication(w, powers::power_of_five_128[index + 1]); 40   4995 value128 secondproduct = full_multiplication(w, powers::power_of_five_128[index + 1]);
HITCBC 41   4995 firstproduct.low += secondproduct.high; 41   4995 firstproduct.low += secondproduct.high;
HITCBC 42   4995 if(secondproduct.high > firstproduct.low) { 42   4995 if(secondproduct.high > firstproduct.low) {
HITCBC 43   1825 firstproduct.high++; 43   1825 firstproduct.high++;
44   } 44   }
45   } 45   }
HITCBC 46   2005879 return firstproduct; 46   2005879 return firstproduct;
47   } 47   }
48   48  
49   namespace detail { 49   namespace detail {
50   /** 50   /**
51   * For q in (0,350), we have that 51   * For q in (0,350), we have that
52   * f = (((152170 + 65536) * q ) >> 16); 52   * f = (((152170 + 65536) * q ) >> 16);
53   * is equal to 53   * is equal to
54   * floor(p) + q 54   * floor(p) + q
55   * where 55   * where
56   * p = log(5**q)/log(2) = q * log(5)/log(2) 56   * p = log(5**q)/log(2) = q * log(5)/log(2)
57   * 57   *
58   * For negative values of q in (-400,0), we have that 58   * For negative values of q in (-400,0), we have that
59   * f = (((152170 + 65536) * q ) >> 16); 59   * f = (((152170 + 65536) * q ) >> 16);
60   * is equal to 60   * is equal to
61   * -ceil(p) + q 61   * -ceil(p) + q
62   * where 62   * where
63   * p = log(5**-q)/log(2) = -q * log(5)/log(2) 63   * p = log(5**-q)/log(2) = -q * log(5)/log(2)
64   */ 64   */
65   constexpr BOOST_FORCEINLINE int32_t power(int32_t q) noexcept { 65   constexpr BOOST_FORCEINLINE int32_t power(int32_t q) noexcept {
HITCBC 66   2005879 return (((152170 + 65536) * q) >> 16) + 63; 66   2005879 return (((152170 + 65536) * q) >> 16) + 63;
67   } 67   }
68   } // namespace detail 68   } // namespace detail
69   69  
70   // create an adjusted mantissa, biased by the invalid power2 70   // create an adjusted mantissa, biased by the invalid power2
71   // for significant digits already multiplied by 10 ** q. 71   // for significant digits already multiplied by 10 ** q.
72   template <typename binary> 72   template <typename binary>
73   BOOST_FORCEINLINE BOOST_JSON_CXX14_CONSTEXPR_NO_INLINE 73   BOOST_FORCEINLINE BOOST_JSON_CXX14_CONSTEXPR_NO_INLINE
74   adjusted_mantissa compute_error_scaled(int64_t q, uint64_t w, int lz) noexcept { 74   adjusted_mantissa compute_error_scaled(int64_t q, uint64_t w, int lz) noexcept {
HITCBC 75   2986 int hilz = int(w >> 63) ^ 1; 75   2986 int hilz = int(w >> 63) ^ 1;
HITCBC 76   2986 adjusted_mantissa answer; 76   2986 adjusted_mantissa answer;
HITCBC 77   2986 answer.mantissa = w << hilz; 77   2986 answer.mantissa = w << hilz;
HITCBC 78   2986 int bias = binary::mantissa_explicit_bits() - binary::minimum_exponent(); 78   2986 int bias = binary::mantissa_explicit_bits() - binary::minimum_exponent();
HITCBC 79   5972 answer.power2 = int32_t(detail::power(int32_t(q)) + bias - hilz - lz - 62 + invalid_am_bias); 79   5972 answer.power2 = int32_t(detail::power(int32_t(q)) + bias - hilz - lz - 62 + invalid_am_bias);
HITCBC 80   2986 return answer; 80   2986 return answer;
81   } 81   }
82   82  
83   // w * 10 ** q, without rounding the representation up. 83   // w * 10 ** q, without rounding the representation up.
84   // the power2 in the exponent will be adjusted by invalid_am_bias. 84   // the power2 in the exponent will be adjusted by invalid_am_bias.
85   template <typename binary> 85   template <typename binary>
86   BOOST_FORCEINLINE BOOST_JSON_FASTFLOAT_CONSTEXPR20 86   BOOST_FORCEINLINE BOOST_JSON_FASTFLOAT_CONSTEXPR20
87   adjusted_mantissa compute_error(int64_t q, uint64_t w) noexcept { 87   adjusted_mantissa compute_error(int64_t q, uint64_t w) noexcept {
HITCBC 88   2986 int lz = leading_zeroes(w); 88   2986 int lz = leading_zeroes(w);
HITCBC 89   2986 w <<= lz; 89   2986 w <<= lz;
90   value128 product = compute_product_approximation<binary::mantissa_explicit_bits() + 3>(q, w); 90   value128 product = compute_product_approximation<binary::mantissa_explicit_bits() + 3>(q, w);
HITCBC 91   8958 return compute_error_scaled<binary>(q, product.high, lz); 91   8958 return compute_error_scaled<binary>(q, product.high, lz);
92   } 92   }
93   93  
94   // w * 10 ** q 94   // w * 10 ** q
95   // The returned value should be a valid ieee64 number that simply need to be packed. 95   // The returned value should be a valid ieee64 number that simply need to be packed.
96   // However, in some very rare cases, the computation will fail. In such cases, we 96   // However, in some very rare cases, the computation will fail. In such cases, we
97   // return an adjusted_mantissa with a negative power of 2: the caller should recompute 97   // return an adjusted_mantissa with a negative power of 2: the caller should recompute
98   // in such cases. 98   // in such cases.
99   template <typename binary> 99   template <typename binary>
100   BOOST_FORCEINLINE BOOST_JSON_FASTFLOAT_CONSTEXPR20 100   BOOST_FORCEINLINE BOOST_JSON_FASTFLOAT_CONSTEXPR20
101   adjusted_mantissa compute_float(int64_t q, uint64_t w) noexcept { 101   adjusted_mantissa compute_float(int64_t q, uint64_t w) noexcept {
HITCBC 102   2004452 adjusted_mantissa answer; 102   2004452 adjusted_mantissa answer;
HITCBC 103   2004452 if ((w == 0) || (q < binary::smallest_power_of_ten())) { 103   2004452 if ((w == 0) || (q < binary::smallest_power_of_ten())) {
HITCBC 104   1 answer.power2 = 0; 104   1 answer.power2 = 0;
HITCBC 105   1 answer.mantissa = 0; 105   1 answer.mantissa = 0;
106   // result should be zero 106   // result should be zero
HITCBC 107   1 return answer; 107   1 return answer;
108   } 108   }
HITCBC 109   2004451 if (q > binary::largest_power_of_ten()) { 109   2004451 if (q > binary::largest_power_of_ten()) {
110   // we want to get infinity: 110   // we want to get infinity:
HITCBC 111   1558 answer.power2 = binary::infinite_power(); 111   1558 answer.power2 = binary::infinite_power();
HITCBC 112   1558 answer.mantissa = 0; 112   1558 answer.mantissa = 0;
HITCBC 113   1558 return answer; 113   1558 return answer;
114   } 114   }
115   // At this point in time q is in [powers::smallest_power_of_five, powers::largest_power_of_five]. 115   // At this point in time q is in [powers::smallest_power_of_five, powers::largest_power_of_five].
116   116  
117   // We want the most significant bit of i to be 1. Shift if needed. 117   // We want the most significant bit of i to be 1. Shift if needed.
HITCBC 118   2002893 int lz = leading_zeroes(w); 118   2002893 int lz = leading_zeroes(w);
HITCBC 119   2002893 w <<= lz; 119   2002893 w <<= lz;
120   120  
121   // The required precision is binary::mantissa_explicit_bits() + 3 because 121   // The required precision is binary::mantissa_explicit_bits() + 3 because
122   // 1. We need the implicit bit 122   // 1. We need the implicit bit
123   // 2. We need an extra bit for rounding purposes 123   // 2. We need an extra bit for rounding purposes
124   // 3. We might lose a bit due to the "upperbit" routine (result too small, requiring a shift) 124   // 3. We might lose a bit due to the "upperbit" routine (result too small, requiring a shift)
125   125  
126   value128 product = compute_product_approximation<binary::mantissa_explicit_bits() + 3>(q, w); 126   value128 product = compute_product_approximation<binary::mantissa_explicit_bits() + 3>(q, w);
127   // The computed 'product' is always sufficient. 127   // The computed 'product' is always sufficient.
128   // Mathematical proof: 128   // Mathematical proof:
129   // Noble Mushtak and Daniel Lemire, Fast Number Parsing Without Fallback (to appear) 129   // Noble Mushtak and Daniel Lemire, Fast Number Parsing Without Fallback (to appear)
130   // See script/mushtak_lemire.py 130   // See script/mushtak_lemire.py
131   131  
132   // The "compute_product_approximation" function can be slightly slower than a branchless approach: 132   // The "compute_product_approximation" function can be slightly slower than a branchless approach:
133   // value128 product = compute_product(q, w); 133   // value128 product = compute_product(q, w);
134   // but in practice, we can win big with the compute_product_approximation if its additional branch 134   // but in practice, we can win big with the compute_product_approximation if its additional branch
135   // is easily predicted. Which is best is data specific. 135   // is easily predicted. Which is best is data specific.
HITCBC 136   2002893 int upperbit = int(product.high >> 63); 136   2002893 int upperbit = int(product.high >> 63);
137   137  
HITCBC 138   2002893 answer.mantissa = product.high >> (upperbit + 64 - binary::mantissa_explicit_bits() - 3); 138   2002893 answer.mantissa = product.high >> (upperbit + 64 - binary::mantissa_explicit_bits() - 3);
139   139  
HITCBC 140   4005786 answer.power2 = int32_t(detail::power(int32_t(q)) + upperbit - lz - binary::minimum_exponent()); 140   4005786 answer.power2 = int32_t(detail::power(int32_t(q)) + upperbit - lz - binary::minimum_exponent());
HITCBC 141   2002893 if (answer.power2 <= 0) { // we have a subnormal? 141   2002893 if (answer.power2 <= 0) { // we have a subnormal?
142   // Here have that answer.power2 <= 0 so -answer.power2 >= 0 142   // Here have that answer.power2 <= 0 so -answer.power2 >= 0
HITCBC 143   82 if(-answer.power2 + 1 >= 64) { // if we have more than 64 bits below the minimum exponent, you have a zero for sure. 143   82 if(-answer.power2 + 1 >= 64) { // if we have more than 64 bits below the minimum exponent, you have a zero for sure.
MISUBC 144   answer.power2 = 0; 144   answer.power2 = 0;
MISUBC 145   answer.mantissa = 0; 145   answer.mantissa = 0;
146   // result should be zero 146   // result should be zero
MISUBC 147   return answer; 147   return answer;
148   } 148   }
149   // next line is safe because -answer.power2 + 1 < 64 149   // next line is safe because -answer.power2 + 1 < 64
HITCBC 150   82 answer.mantissa >>= -answer.power2 + 1; 150   82 answer.mantissa >>= -answer.power2 + 1;
151   // Thankfully, we can't have both "round-to-even" and subnormals because 151   // Thankfully, we can't have both "round-to-even" and subnormals because
152   // "round-to-even" only occurs for powers close to 0. 152   // "round-to-even" only occurs for powers close to 0.
HITCBC 153   82 answer.mantissa += (answer.mantissa & 1); // round up 153   82 answer.mantissa += (answer.mantissa & 1); // round up
HITCBC 154   82 answer.mantissa >>= 1; 154   82 answer.mantissa >>= 1;
155   // There is a weird scenario where we don't have a subnormal but just. 155   // There is a weird scenario where we don't have a subnormal but just.
156   // Suppose we start with 2.2250738585072013e-308, we end up 156   // Suppose we start with 2.2250738585072013e-308, we end up
157   // with 0x3fffffffffffff x 2^-1023-53 which is technically subnormal 157   // with 0x3fffffffffffff x 2^-1023-53 which is technically subnormal
158   // whereas 0x40000000000000 x 2^-1023-53 is normal. Now, we need to round 158   // whereas 0x40000000000000 x 2^-1023-53 is normal. Now, we need to round
159   // up 0x3fffffffffffff x 2^-1023-53 and once we do, we are no longer 159   // up 0x3fffffffffffff x 2^-1023-53 and once we do, we are no longer
160   // subnormal, but we can only know this after rounding. 160   // subnormal, but we can only know this after rounding.
161   // So we only declare a subnormal if we are smaller than the threshold. 161   // So we only declare a subnormal if we are smaller than the threshold.
HITCBC 162   82 answer.power2 = (answer.mantissa < (uint64_t(1) << binary::mantissa_explicit_bits())) ? 0 : 1; 162   82 answer.power2 = (answer.mantissa < (uint64_t(1) << binary::mantissa_explicit_bits())) ? 0 : 1;
HITCBC 163   82 return answer; 163   82 return answer;
164   } 164   }
165   165  
166   // usually, we round *up*, but if we fall right in between and and we have an 166   // usually, we round *up*, but if we fall right in between and and we have an
167   // even basis, we need to round down 167   // even basis, we need to round down
168   // We are only concerned with the cases where 5**q fits in single 64-bit word. 168   // We are only concerned with the cases where 5**q fits in single 64-bit word.
HITCBC 169   2009538 if ((product.low <= 1) && (q >= binary::min_exponent_round_to_even()) && (q <= binary::max_exponent_round_to_even()) && 169   2009538 if ((product.low <= 1) && (q >= binary::min_exponent_round_to_even()) && (q <= binary::max_exponent_round_to_even()) &&
HITCBC 170   6727 ((answer.mantissa & 3) == 1) ) { // we may fall between two floats! 170   6727 ((answer.mantissa & 3) == 1) ) { // we may fall between two floats!
171   // To be in-between two floats we need that in doing 171   // To be in-between two floats we need that in doing
172   // answer.mantissa = product.high >> (upperbit + 64 - binary::mantissa_explicit_bits() - 3); 172   // answer.mantissa = product.high >> (upperbit + 64 - binary::mantissa_explicit_bits() - 3);
173   // ... we dropped out only zeroes. But if this happened, then we can go back!!! 173   // ... we dropped out only zeroes. But if this happened, then we can go back!!!
HITCBC 174   1857 if((answer.mantissa << (upperbit + 64 - binary::mantissa_explicit_bits() - 3)) == product.high) { 174   1857 if((answer.mantissa << (upperbit + 64 - binary::mantissa_explicit_bits() - 3)) == product.high) {
HITCBC 175   14 answer.mantissa &= ~uint64_t(1); // flip it so that we do not round up 175   14 answer.mantissa &= ~uint64_t(1); // flip it so that we do not round up
176   } 176   }
177   } 177   }
178   178  
HITCBC 179   2002811 answer.mantissa += (answer.mantissa & 1); // round up 179   2002811 answer.mantissa += (answer.mantissa & 1); // round up
HITCBC 180   2002811 answer.mantissa >>= 1; 180   2002811 answer.mantissa >>= 1;
HITCBC 181   2002811 if (answer.mantissa >= (uint64_t(2) << binary::mantissa_explicit_bits())) { 181   2002811 if (answer.mantissa >= (uint64_t(2) << binary::mantissa_explicit_bits())) {
HITCBC 182   360 answer.mantissa = (uint64_t(1) << binary::mantissa_explicit_bits()); 182   360 answer.mantissa = (uint64_t(1) << binary::mantissa_explicit_bits());
HITCBC 183   360 answer.power2++; // undo previous addition 183   360 answer.power2++; // undo previous addition
184   } 184   }
185   185  
HITCBC 186   2002811 answer.mantissa &= ~(uint64_t(1) << binary::mantissa_explicit_bits()); 186   2002811 answer.mantissa &= ~(uint64_t(1) << binary::mantissa_explicit_bits());
HITCBC 187   2002811 if (answer.power2 >= binary::infinite_power()) { // infinity 187   2002811 if (answer.power2 >= binary::infinite_power()) { // infinity
HITCBC 188   59658 answer.power2 = binary::infinite_power(); 188   59658 answer.power2 = binary::infinite_power();
HITCBC 189   59658 answer.mantissa = 0; 189   59658 answer.mantissa = 0;
190   } 190   }
HITCBC 191   2002811 return answer; 191   2002811 return answer;
192   } 192   }
193   193  
194   }}}}}} // namespace fast_float 194   }}}}}} // namespace fast_float
195   195  
196   #endif 196   #endif