|
| 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 |
0 commit comments