| Line | Branch | Exec | Source |
|---|---|---|---|
| 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 | #ifndef BOOST_JSON_DETAIL_CHARCONV_DETAIL_COMPUTE_FLOAT64_HPP | ||
| 7 | #define BOOST_JSON_DETAIL_CHARCONV_DETAIL_COMPUTE_FLOAT64_HPP | ||
| 8 | |||
| 9 | #include <boost/json/detail/charconv/detail/config.hpp> | ||
| 10 | #include <boost/json/detail/charconv/detail/significand_tables.hpp> | ||
| 11 | #include <boost/json/detail/charconv/detail/emulated128.hpp> | ||
| 12 | #include <boost/core/bit.hpp> | ||
| 13 | #include <cstdint> | ||
| 14 | #include <cfloat> | ||
| 15 | #include <cstring> | ||
| 16 | #include <cmath> | ||
| 17 | |||
| 18 | namespace boost { namespace json { namespace detail { namespace charconv { namespace detail { | ||
| 19 | |||
| 20 | static constexpr double powers_of_ten[] = { | ||
| 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 | ||
| 23 | }; | ||
| 24 | |||
| 25 | // Attempts to compute i * 10^(power) exactly; and if "negative" is true, negate the result. | ||
| 26 | // | ||
| 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). | ||
| 29 | // We assume that power is in the [-325, 308] interval. | ||
| 30 | ✗ | inline double compute_float64(std::int64_t power, std::uint64_t i, bool negative, bool& success) noexcept | |
| 31 | { | ||
| 32 | static constexpr auto smallest_power = -325; | ||
| 33 | static constexpr auto largest_power = 308; | ||
| 34 | |||
| 35 | // We start with a fast path | ||
| 36 | // It was described in Clinger WD. | ||
| 37 | // How to read floating point numbers accurately. | ||
| 38 | // ACM SIGPLAN Notices. 1990 | ||
| 39 | #if (FLT_EVAL_METHOD != 1) && (FLT_EVAL_METHOD != 0) | ||
| 40 | if (0 <= power && power <= 22 && i <= UINT64_C(9007199254740991)) | ||
| 41 | #else | ||
| 42 | ✗ | if (-22 <= power && power <= 22 && i <= UINT64_C(9007199254740991)) | |
| 43 | #endif | ||
| 44 | { | ||
| 45 | // The general idea is as follows. | ||
| 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 | ||
| 48 | // values | ||
| 49 | // (binary64). | ||
| 50 | // 2) Because s and p can be represented exactly as floating-point values, | ||
| 51 | // then s * p | ||
| 52 | // and s / p will produce correctly rounded values. | ||
| 53 | |||
| 54 | ✗ | auto d = static_cast<double>(i); | |
| 55 | |||
| 56 | ✗ | if (power < 0) | |
| 57 | { | ||
| 58 | ✗ | d = d / powers_of_ten[-power]; | |
| 59 | } | ||
| 60 | else | ||
| 61 | { | ||
| 62 | ✗ | d = d * powers_of_ten[power]; | |
| 63 | } | ||
| 64 | |||
| 65 | ✗ | if (negative) | |
| 66 | { | ||
| 67 | ✗ | d = -d; | |
| 68 | } | ||
| 69 | |||
| 70 | ✗ | success = true; | |
| 71 | ✗ | return d; | |
| 72 | } | ||
| 73 | |||
| 74 | // When 22 < power && power < 22 + 16, we could | ||
| 75 | // hope for another, secondary fast path. It was | ||
| 76 | // described by David M. Gay in "Correctly rounded | ||
| 77 | // binary-decimal and decimal-binary conversions." (1990) | ||
| 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 | ||
| 80 | // (e.g., when i * 10^x < 2^53), | ||
| 81 | // then you can still proceed and do (i * 10^x) * 10^22. | ||
| 82 | // Is this worth your time? | ||
| 83 | // You need 22 < power *and* power < 22 + 16 *and* (i * 10^(x-22) < 2^53) | ||
| 84 | // for this second fast path to work. | ||
| 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 | ||
| 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: | ||
| 89 | // http://www.exploringbinary.com/fast-path-decimal-to-floating-point-conversion/ | ||
| 90 | // also used in RapidJSON: https://rapidjson.org/strtod_8h_source.html | ||
| 91 | |||
| 92 | ✗ | if (i == 0 || power < smallest_power) | |
| 93 | { | ||
| 94 | ✗ | return negative ? -0.0 : 0.0; | |
| 95 | } | ||
| 96 | ✗ | else if (power > largest_power) | |
| 97 | { | ||
| 98 | ✗ | return negative ? -HUGE_VAL : HUGE_VAL; | |
| 99 | } | ||
| 100 | |||
| 101 | ✗ | const std::uint64_t factor_significand = significand_64[power - smallest_power]; | |
| 102 | ✗ | const std::int64_t exponent = (((152170 + 65536) * power) >> 16) + 1024 + 63; | |
| 103 | ✗ | int leading_zeros = boost::core::countl_zero(i); | |
| 104 | ✗ | i <<= static_cast<std::uint64_t>(leading_zeros); | |
| 105 | |||
| 106 | ✗ | uint128 product = umul128(i, factor_significand); | |
| 107 | ✗ | std::uint64_t low = product.low; | |
| 108 | ✗ | std::uint64_t high = product.high; | |
| 109 | |||
| 110 | // We know that upper has at most one leading zero because | ||
| 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). | ||
| 113 | // | ||
| 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 | ||
| 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, | ||
| 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. | ||
| 120 | // | ||
| 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 | ||
| 123 | // lower + i < lower to be true (proba. much higher than 1%). | ||
| 124 | ✗ | if (BOOST_UNLIKELY((high & 0x1FF) == 0x1FF) && (low + i < low)) | |
| 125 | { | ||
| 126 | ✗ | const std::uint64_t factor_significand_low = significand_128[power - smallest_power]; | |
| 127 | ✗ | product = umul128(i, factor_significand_low); | |
| 128 | //const std::uint64_t product_low = product.low; | ||
| 129 | ✗ | const std::uint64_t product_middle2 = product.high; | |
| 130 | ✗ | const std::uint64_t product_middle1 = low; | |
| 131 | ✗ | std::uint64_t product_high = high; | |
| 132 | ✗ | const std::uint64_t product_middle = product_middle1 + product_middle2; | |
| 133 | |||
| 134 | ✗ | if (product_middle < product_middle1) | |
| 135 | { | ||
| 136 | ✗ | product_high++; | |
| 137 | } | ||
| 138 | |||
| 139 | // Commented out because possibly unneeded | ||
| 140 | // See: https://arxiv.org/pdf/2212.06644.pdf | ||
| 141 | /* | ||
| 142 | // we want to check whether mantissa *i + i would affect our result | ||
| 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))) | ||
| 145 | { | ||
| 146 | success = false; | ||
| 147 | return 0; | ||
| 148 | } | ||
| 149 | */ | ||
| 150 | |||
| 151 | ✗ | low = product_middle; | |
| 152 | ✗ | high = product_high; | |
| 153 | } | ||
| 154 | |||
| 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 | ||
| 157 | ✗ | const std::uint64_t upper_bit = high >> 63; | |
| 158 | ✗ | std::uint64_t significand = high >> (upper_bit + 9); | |
| 159 | ✗ | leading_zeros += static_cast<int>(1 ^ upper_bit); | |
| 160 | |||
| 161 | // If we have lots of trailing zeros we may fall between two values | ||
| 162 | ✗ | if (BOOST_UNLIKELY((low == 0) && ((high & 0x1FF) == 0) && ((significand & 3) == 1))) | |
| 163 | { | ||
| 164 | // if significand & 1 == 1 we might need to round up | ||
| 165 | ✗ | success = false; | |
| 166 | ✗ | return 0; | |
| 167 | } | ||
| 168 | |||
| 169 | ✗ | significand += significand & 1; | |
| 170 | ✗ | significand >>= 1; | |
| 171 | |||
| 172 | // Here the significand < (1<<53), unless there is an overflow | ||
| 173 | ✗ | if (significand >= (UINT64_C(1) << 53)) | |
| 174 | { | ||
| 175 | ✗ | significand = (UINT64_C(1) << 52); | |
| 176 | ✗ | leading_zeros--; | |
| 177 | } | ||
| 178 | |||
| 179 | ✗ | significand &= ~(UINT64_C(1) << 52); | |
| 180 | ✗ | const std::uint64_t real_exponent = exponent - leading_zeros; | |
| 181 | |||
| 182 | // We have to check that real_exponent is in range, otherwise fail | ||
| 183 | ✗ | if (BOOST_UNLIKELY((real_exponent < 1) || (real_exponent > 2046))) | |
| 184 | { | ||
| 185 | ✗ | success = false; | |
| 186 | ✗ | return 0; | |
| 187 | } | ||
| 188 | |||
| 189 | ✗ | significand |= real_exponent << 52; | |
| 190 | ✗ | significand |= ((static_cast<std::uint64_t>(negative) << 63)); | |
| 191 | |||
| 192 | double d; | ||
| 193 | ✗ | std::memcpy(&d, &significand, sizeof(d)); | |
| 194 | |||
| 195 | ✗ | success = true; | |
| 196 | ✗ | return d; | |
| 197 | } | ||
| 198 | |||
| 199 | }}}}} // Namespaces | ||
| 200 | |||
| 201 | #endif // BOOST_JSON_DETAIL_CHARCONV_DETAIL_COMPUTE_FLOAT64_HPP | ||
| 202 |