include/boost/json /detail /charconv /detail /compute_float64.hpp
0.00% Lines (0/53) 0.00% Functions (0/1)
[ ] prev/next file n p prev/next uncovered line
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   #ifndef BOOST_JSON_DETAIL_CHARCONV_DETAIL_COMPUTE_FLOAT64_HPP 6   #ifndef BOOST_JSON_DETAIL_CHARCONV_DETAIL_COMPUTE_FLOAT64_HPP
7   #define BOOST_JSON_DETAIL_CHARCONV_DETAIL_COMPUTE_FLOAT64_HPP 7   #define BOOST_JSON_DETAIL_CHARCONV_DETAIL_COMPUTE_FLOAT64_HPP
8   8  
9   #include <boost/json/detail/charconv/detail/config.hpp> 9   #include <boost/json/detail/charconv/detail/config.hpp>
10   #include <boost/json/detail/charconv/detail/significand_tables.hpp> 10   #include <boost/json/detail/charconv/detail/significand_tables.hpp>
11   #include <boost/json/detail/charconv/detail/emulated128.hpp> 11   #include <boost/json/detail/charconv/detail/emulated128.hpp>
12   #include <boost/core/bit.hpp> 12   #include <boost/core/bit.hpp>
13   #include <cstdint> 13   #include <cstdint>
14   #include <cfloat> 14   #include <cfloat>
15   #include <cstring> 15   #include <cstring>
16   #include <cmath> 16   #include <cmath>
17   17  
18   namespace boost { namespace json { namespace detail { namespace charconv { namespace detail { 18   namespace boost { namespace json { namespace detail { namespace charconv { namespace detail {
19   19  
20   static constexpr double powers_of_ten[] = { 20   static constexpr double powers_of_ten[] = {
21   1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9, 1e10, 1e11, 21   1e0, 1e1, 1e2, 1e3, 1e4, 1e5, 1e6, 1e7, 1e8, 1e9, 1e10, 1e11,
22   1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19, 1e20, 1e21, 1e22 22   1e12, 1e13, 1e14, 1e15, 1e16, 1e17, 1e18, 1e19, 1e20, 1e21, 1e22
23   }; 23   };
24   24  
25   // Attempts to compute i * 10^(power) exactly; and if "negative" is true, negate the result. 25   // Attempts to compute i * 10^(power) exactly; and if "negative" is true, negate the result.
26   // 26   //
27   // This function will only work in some cases, when it does not work, success is 27   // This function will only work in some cases, when it does not work, success is
28   // set to false. This should work *most of the time* (like 99% of the time). 28   // set to false. This should work *most of the time* (like 99% of the time).
29   // We assume that power is in the [-325, 308] interval. 29   // We assume that power is in the [-325, 308] interval.
MISUBC 30   ✗ inline double compute_float64(std::int64_t power, std::uint64_t i, bool negative, bool& success) noexcept 30   ✗ inline double compute_float64(std::int64_t power, std::uint64_t i, bool negative, bool& success) noexcept
31   { 31   {
32   static constexpr auto smallest_power = -325; 32   static constexpr auto smallest_power = -325;
33   static constexpr auto largest_power = 308; 33   static constexpr auto largest_power = 308;
34   34  
35   // We start with a fast path 35   // We start with a fast path
36   // It was described in Clinger WD. 36   // It was described in Clinger WD.
37   // How to read floating point numbers accurately. 37   // How to read floating point numbers accurately.
38   // ACM SIGPLAN Notices. 1990 38   // ACM SIGPLAN Notices. 1990
39   #if (FLT_EVAL_METHOD != 1) && (FLT_EVAL_METHOD != 0) 39   #if (FLT_EVAL_METHOD != 1) && (FLT_EVAL_METHOD != 0)
40   if (0 <= power && power <= 22 && i <= UINT64_C(9007199254740991)) 40   if (0 <= power && power <= 22 && i <= UINT64_C(9007199254740991))
41   #else 41   #else
MISUBC 42   ✗ if (-22 <= power && power <= 22 && i <= UINT64_C(9007199254740991)) 42   ✗ if (-22 <= power && power <= 22 && i <= UINT64_C(9007199254740991))
43   #endif 43   #endif
44   { 44   {
45   // The general idea is as follows. 45   // The general idea is as follows.
46   // If 0 <= s < 2^53 and if 10^0 <= p <= 10^22 then 46   // If 0 <= s < 2^53 and if 10^0 <= p <= 10^22 then
47   // 1) Both s and p can be represented exactly as 64-bit floating-point 47   // 1) Both s and p can be represented exactly as 64-bit floating-point
48   // values 48   // values
49   // (binary64). 49   // (binary64).
50   // 2) Because s and p can be represented exactly as floating-point values, 50   // 2) Because s and p can be represented exactly as floating-point values,
51   // then s * p 51   // then s * p
52   // and s / p will produce correctly rounded values. 52   // and s / p will produce correctly rounded values.
53   53  
MISUBC 54   ✗ auto d = static_cast<double>(i); 54   ✗ auto d = static_cast<double>(i);
55   55  
MISUBC 56   ✗ if (power < 0) 56   ✗ if (power < 0)
57   { 57   {
MISUBC 58   ✗ d = d / powers_of_ten[-power]; 58   ✗ d = d / powers_of_ten[-power];
59   } 59   }
60   else 60   else
61   { 61   {
MISUBC 62   ✗ d = d * powers_of_ten[power]; 62   ✗ d = d * powers_of_ten[power];
63   } 63   }
64   64  
MISUBC 65   ✗ if (negative) 65   ✗ if (negative)
66   { 66   {
MISUBC 67   ✗ d = -d; 67   ✗ d = -d;
68   } 68   }
69   69  
MISUBC 70   ✗ success = true; 70   ✗ success = true;
MISUBC 71   ✗ return d; 71   ✗ return d;
72   } 72   }
73   73  
74   // When 22 < power && power < 22 + 16, we could 74   // When 22 < power && power < 22 + 16, we could
75   // hope for another, secondary fast path. It was 75   // hope for another, secondary fast path. It was
76   // described by David M. Gay in "Correctly rounded 76   // described by David M. Gay in "Correctly rounded
77   // binary-decimal and decimal-binary conversions." (1990) 77   // binary-decimal and decimal-binary conversions." (1990)
78   // If you need to compute i * 10^(22 + x) for x < 16, 78   // If you need to compute i * 10^(22 + x) for x < 16,
79   // first compute i * 10^x, if you know that result is exact 79   // first compute i * 10^x, if you know that result is exact
80   // (e.g., when i * 10^x < 2^53), 80   // (e.g., when i * 10^x < 2^53),
81   // then you can still proceed and do (i * 10^x) * 10^22. 81   // then you can still proceed and do (i * 10^x) * 10^22.
82   // Is this worth your time? 82   // Is this worth your time?
83   // You need 22 < power *and* power < 22 + 16 *and* (i * 10^(x-22) < 2^53) 83   // You need 22 < power *and* power < 22 + 16 *and* (i * 10^(x-22) < 2^53)
84   // for this second fast path to work. 84   // for this second fast path to work.
85   // If you have 22 < power *and* power < 22 + 16, and then you 85   // If you have 22 < power *and* power < 22 + 16, and then you
86   // optimistically compute "i * 10^(x-22)", there is still a chance that you 86   // optimistically compute "i * 10^(x-22)", there is still a chance that you
87   // have wasted your time if i * 10^(x-22) >= 2^53. It makes the use cases of 87   // have wasted your time if i * 10^(x-22) >= 2^53. It makes the use cases of
88   // this optimization maybe less common than we would like. Source: 88   // this optimization maybe less common than we would like. Source:
89   // http://www.exploringbinary.com/fast-path-decimal-to-floating-point-conversion/ 89   // http://www.exploringbinary.com/fast-path-decimal-to-floating-point-conversion/
90   // also used in RapidJSON: https://rapidjson.org/strtod_8h_source.html 90   // also used in RapidJSON: https://rapidjson.org/strtod_8h_source.html
91   91  
MISUBC 92   ✗ if (i == 0 || power < smallest_power) 92   ✗ if (i == 0 || power < smallest_power)
93   { 93   {
MISUBC 94   ✗ return negative ? -0.0 : 0.0; 94   ✗ return negative ? -0.0 : 0.0;
95   } 95   }
MISUBC 96   ✗ else if (power > largest_power) 96   ✗ else if (power > largest_power)
97   { 97   {
MISUBC 98   ✗ return negative ? -HUGE_VAL : HUGE_VAL; 98   ✗ return negative ? -HUGE_VAL : HUGE_VAL;
99   } 99   }
100   100  
MISUBC 101   ✗ const std::uint64_t factor_significand = significand_64[power - smallest_power]; 101   ✗ const std::uint64_t factor_significand = significand_64[power - smallest_power];
MISUBC 102   ✗ const std::int64_t exponent = (((152170 + 65536) * power) >> 16) + 1024 + 63; 102   ✗ const std::int64_t exponent = (((152170 + 65536) * power) >> 16) + 1024 + 63;
MISUBC 103   ✗ int leading_zeros = boost::core::countl_zero(i); 103   ✗ int leading_zeros = boost::core::countl_zero(i);
MISUBC 104   ✗ i <<= static_cast<std::uint64_t>(leading_zeros); 104   ✗ i <<= static_cast<std::uint64_t>(leading_zeros);
105   105  
MISUBC 106   ✗ uint128 product = umul128(i, factor_significand); 106   ✗ uint128 product = umul128(i, factor_significand);
MISUBC 107   ✗ std::uint64_t low = product.low; 107   ✗ std::uint64_t low = product.low;
MISUBC 108   ✗ std::uint64_t high = product.high; 108   ✗ std::uint64_t high = product.high;
109   109  
110   // We know that upper has at most one leading zero because 110   // We know that upper has at most one leading zero because
111   // both i and factor_mantissa have a leading one. This means 111   // both i and factor_mantissa have a leading one. This means
112   // that the result is at least as large as ((1<<63)*(1<<63))/(1<<64). 112   // that the result is at least as large as ((1<<63)*(1<<63))/(1<<64).
113   // 113   //
114   // As long as the first 9 bits of "upper" are not "1", then we 114   // As long as the first 9 bits of "upper" are not "1", then we
115   // know that we have an exact computed value for the leading 115   // know that we have an exact computed value for the leading
116   // 55 bits because any imprecision would play out as a +1, in the worst case. 116   // 55 bits because any imprecision would play out as a +1, in the worst case.
117   // Having 55 bits is necessary because we need 53 bits for the mantissa, 117   // Having 55 bits is necessary because we need 53 bits for the mantissa,
118   // but we have to have one rounding bit and, we can waste a bit if the most 118   // but we have to have one rounding bit and, we can waste a bit if the most
119   // significant bit of the product is zero. 119   // significant bit of the product is zero.
120   // 120   //
121   // We expect this next branch to be rarely taken (say 1% of the time). 121   // We expect this next branch to be rarely taken (say 1% of the time).
122   // When (upper & 0x1FF) == 0x1FF, it can be common for 122   // When (upper & 0x1FF) == 0x1FF, it can be common for
123   // lower + i < lower to be true (proba. much higher than 1%). 123   // lower + i < lower to be true (proba. much higher than 1%).
MISUBC 124   ✗ if (BOOST_UNLIKELY((high & 0x1FF) == 0x1FF) && (low + i < low)) 124   ✗ if (BOOST_UNLIKELY((high & 0x1FF) == 0x1FF) && (low + i < low))
125   { 125   {
MISUBC 126   ✗ const std::uint64_t factor_significand_low = significand_128[power - smallest_power]; 126   ✗ const std::uint64_t factor_significand_low = significand_128[power - smallest_power];
MISUBC 127   ✗ product = umul128(i, factor_significand_low); 127   ✗ product = umul128(i, factor_significand_low);
128   //const std::uint64_t product_low = product.low; 128   //const std::uint64_t product_low = product.low;
MISUBC 129   ✗ const std::uint64_t product_middle2 = product.high; 129   ✗ const std::uint64_t product_middle2 = product.high;
MISUBC 130   ✗ const std::uint64_t product_middle1 = low; 130   ✗ const std::uint64_t product_middle1 = low;
MISUBC 131   ✗ std::uint64_t product_high = high; 131   ✗ std::uint64_t product_high = high;
MISUBC 132   ✗ const std::uint64_t product_middle = product_middle1 + product_middle2; 132   ✗ const std::uint64_t product_middle = product_middle1 + product_middle2;
133   133  
MISUBC 134   ✗ if (product_middle < product_middle1) 134   ✗ if (product_middle < product_middle1)
135   { 135   {
MISUBC 136   ✗ product_high++; 136   ✗ product_high++;
137   } 137   }
138   138  
139   // Commented out because possibly unneeded 139   // Commented out because possibly unneeded
140   // See: https://arxiv.org/pdf/2212.06644.pdf 140   // See: https://arxiv.org/pdf/2212.06644.pdf
141   /* 141   /*
142   // we want to check whether mantissa *i + i would affect our result 142   // we want to check whether mantissa *i + i would affect our result
143   // This does happen, e.g. with 7.3177701707893310e+15 143   // This does happen, e.g. with 7.3177701707893310e+15
144   if (((product_middle + 1 == 0) && ((product_high & 0x1FF) == 0x1FF) && (product_low + i < product_low))) 144   if (((product_middle + 1 == 0) && ((product_high & 0x1FF) == 0x1FF) && (product_low + i < product_low)))
145   { 145   {
146   success = false; 146   success = false;
147   return 0; 147   return 0;
148   } 148   }
149   */ 149   */
150   150  
MISUBC 151   ✗ low = product_middle; 151   ✗ low = product_middle;
MISUBC 152   ✗ high = product_high; 152   ✗ high = product_high;
153   } 153   }
154   154  
155   // The final significand should be 53 bits with a leading 1 155   // The final significand should be 53 bits with a leading 1
156   // We shift it so that it occupies 54 bits with a leading 1 156   // We shift it so that it occupies 54 bits with a leading 1
MISUBC 157   ✗ const std::uint64_t upper_bit = high >> 63; 157   ✗ const std::uint64_t upper_bit = high >> 63;
MISUBC 158   ✗ std::uint64_t significand = high >> (upper_bit + 9); 158   ✗ std::uint64_t significand = high >> (upper_bit + 9);
MISUBC 159   ✗ leading_zeros += static_cast<int>(1 ^ upper_bit); 159   ✗ leading_zeros += static_cast<int>(1 ^ upper_bit);
160   160  
161   // If we have lots of trailing zeros we may fall between two values 161   // If we have lots of trailing zeros we may fall between two values
MISUBC 162   ✗ if (BOOST_UNLIKELY((low == 0) && ((high & 0x1FF) == 0) && ((significand & 3) == 1))) 162   ✗ if (BOOST_UNLIKELY((low == 0) && ((high & 0x1FF) == 0) && ((significand & 3) == 1)))
163   { 163   {
164   // if significand & 1 == 1 we might need to round up 164   // if significand & 1 == 1 we might need to round up
MISUBC 165   ✗ success = false; 165   ✗ success = false;
MISUBC 166   ✗ return 0; 166   ✗ return 0;
167   } 167   }
168   168  
MISUBC 169   ✗ significand += significand & 1; 169   ✗ significand += significand & 1;
MISUBC 170   ✗ significand >>= 1; 170   ✗ significand >>= 1;
171   171  
172   // Here the significand < (1<<53), unless there is an overflow 172   // Here the significand < (1<<53), unless there is an overflow
MISUBC 173   ✗ if (significand >= (UINT64_C(1) << 53)) 173   ✗ if (significand >= (UINT64_C(1) << 53))
174   { 174   {
MISUBC 175   ✗ significand = (UINT64_C(1) << 52); 175   ✗ significand = (UINT64_C(1) << 52);
MISUBC 176   ✗ leading_zeros--; 176   ✗ leading_zeros--;
177   } 177   }
178   178  
MISUBC 179   ✗ significand &= ~(UINT64_C(1) << 52); 179   ✗ significand &= ~(UINT64_C(1) << 52);
MISUBC 180   ✗ const std::uint64_t real_exponent = exponent - leading_zeros; 180   ✗ const std::uint64_t real_exponent = exponent - leading_zeros;
181   181  
182   // We have to check that real_exponent is in range, otherwise fail 182   // We have to check that real_exponent is in range, otherwise fail
MISUBC 183   ✗ if (BOOST_UNLIKELY((real_exponent < 1) || (real_exponent > 2046))) 183   ✗ if (BOOST_UNLIKELY((real_exponent < 1) || (real_exponent > 2046)))
184   { 184   {
MISUBC 185   ✗ success = false; 185   ✗ success = false;
MISUBC 186   ✗ return 0; 186   ✗ return 0;
187   } 187   }
188   188  
MISUBC 189   ✗ significand |= real_exponent << 52; 189   ✗ significand |= real_exponent << 52;
MISUBC 190   ✗ significand |= ((static_cast<std::uint64_t>(negative) << 63)); 190   ✗ significand |= ((static_cast<std::uint64_t>(negative) << 63));
191   191  
192   double d; 192   double d;
MISUBC 193   ✗ std::memcpy(&d, &significand, sizeof(d)); 193   ✗ std::memcpy(&d, &significand, sizeof(d));
194   194  
MISUBC 195   ✗ success = true; 195   ✗ success = true;
MISUBC 196   ✗ return d; 196   ✗ return d;
197   } 197   }
198   198  
199   }}}}} // Namespaces 199   }}}}} // Namespaces
200   200  
201   #endif // BOOST_JSON_DETAIL_CHARCONV_DETAIL_COMPUTE_FLOAT64_HPP 201   #endif // BOOST_JSON_DETAIL_CHARCONV_DETAIL_COMPUTE_FLOAT64_HPP